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V/ave  propagation  to  far-regional  distances  in  the  Western 

United  States 

Karl  Koch,  Briar\  Stump 
Southern  Methodist  University 
Department  of  Geological  Sciences 
Dallas,  TX  75275 


1.  Introduction 

Recent  work  in  the  field  of  discrimination,  verification,  and  yield  estimation  has 
been  directed  towards  utilization  of  regional  seismograms  for  a  number  of 
reasons.  First,  these  data  may  be  considered  as  the  basis  for  verification  in  the 
context  of  a  reduced  threshold  test  ban  treaty.  Second,  one  regional  phase,  Lg,  has 
proven  to  be  not  only  a  stable  measure  for  yield  (Nuttli  1986;  Patton  1988),  but  is 
itself  a  useful  waveform  feature,  as  it  is  often  associated  with  the  largest 
amplitudes  observed  in  regional  seismograms.  This  is  particularly  true  for  small 
events,  where  Lg  in  many  cases  (Hansen  et.al.  1990)  is  the  only  phase  that  can  be 
identified  from  seismic  observations. 

Using  regional  phases  for  discrimination  and  yield  estimation  has  been  quite 
successful  in  recent  studies.  Nuttli  (1986a,b;  1988;  Patton  1988)  found  that  Lg  is  a 
very  stable  yield  estimator  and  that  the  mbLg  yield  estimator  can  be  easily 
transported  to  different  geological  settings.  Taylor  &  Randall  (1988)  reported  the 
usefulness  of  spectral  measurement  from  regional  phase  (spectral  ratios)  for 
discrimination  purposes.  Other  authors  (Pomeroy  et  al.  1982;  Bennett  &  Murphy 
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1986)  identify  amplitude  as  well  as  spectral  ratios  among  different  regional 
phases  as  useful  discriminants.  All  these  empirical  results  demand  a  more 
thorough  theoretical  evaluation  of  wave  propagation  to  regional  distances  in 
order  to  provide  a  physical  understanding  for  these  measures  as  well  as  separate 
the  source  contributions  from  the  propagation  path  effects.  The  work  by  Lilwall 
(1988)  addressed  this  problem  by  studying  different  discriminants  with  synthetic 
models  of  regional  seismograms.  Barker  et.al  (1990)  used  a  similar  approach  to 
investigate  the  source  contributions  for  a  variety  of  source  types  on  particular 
regional  phases,  especially  the  relative  excitation  of  these  phases.  Neither  of 
these  authors,  however,  matched  regional  seismograms  to  observed  recordings. 
This  study  however  will  try  to  fill  this  gap  by  attempting  to  attribute  as  much  as 
possible  of  observed  waveforms  characteristics  to  propagation  path  effects. 
Following  this  study  we  will  address  differences  in  regional  phases  that  are 
related  to  different  source  representation  accompanying  the  nuclear  explosion. 

The  distance  range  for  the  regional  data  of  this  study  is  from  800-2000  km  across 
the  Basin  and  Range  province  (Western  United  States),  which  is  effectively  what 
can  be  expected  for  the  station  spacing  in  host  countries  under  a  verification 
scenario.  This  distance  range  implies  that  the  region  illuminated  by  waves 
travelling  from  the  source  to  the  receiver  is  the  crust  extending  into  the  upper 
mantle.  A  number  of  velocity  models  for  crustal  and  upper  mantle  propagation 
have  been  developed  for  the  Western  U.tited  States.  Many  of  these  models  have 
been  derived  by  analysis  of  long-period  body  waves  (Archambeau  et  al  1969, 
Johnson  1967,  Wiggins  &  Helmberger  1973,  Burdick  &  Helmberger  1978, 
Helmberger  &  Engen  1974)  and  from  surface  waves  (Priestley  &  Brune  1978, 
Patton  &  Taylor  1984).  None  of  these  modeling  studies  have  tried  to  incorporate 
complete  waveforms  extending  from  the  body  waves  to  the  surface  waves.  The 
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focus  of  this  study  will  be  to  integrate  all  the  observed  phases  in  the  regional 
seismograms. 

The  characteristics  of  previous  models  include  a  simple  crust,  with  upper  crustal 
compressional  wave  velocities  arotmd  6  km/ s  and  a  lower  crust  with  a  P  velocity 
of  about  6.5  km/s.  The  thickness  of  the  crust  is  35±5  km.  The  upper  mantle  is 
often  associated  with  a  high-velocity  structure  below  the  Moho,  underlain  by  a 
low  velocity  region  (both  in  compressional  and  shear  wave  velocity).  This  low 
velocity  region  extends  for  P  waves  frequently  to  a  depth  of  about  150,  while  for 
shear  waves  the  zone  of  low  S  velocities  extends  to  somewhat  greater  depth. 
Most  P  models  incorporate  a  discontinuous  decrease  of  velocity  into  the  low 
velocity  zone  followed  by  a  large  positive  gradient,  while  the  S  wave  speed  in  the 
low  velocity  zone  is  characterized  by  an  almost  constant  velocity  throughout  the 
region.  As  is  indicated  by  our  data,  the  velocity  model  below  about  250  km  is  not 
important  in  terms  of  wave  propagations  effects  for  the  distance  range 
considered. 

Since  Olsen  et  al,  (1980,  1983)  dealt  with  crustal  and  upper  mantle  wave 
propagation  to  distances  between  600  and  900  km,  we  use  these  velocity  models 
as  the  starting  point  in  synthesizing  the  characteristics  of  observed  waveforms 
which  will  be  discussed  in  the  following  section..  The  analysis  of  the  far-regional 
waveforms  will  focus  on  the  major  regional  phases,  which  include  body  wave 
arrivals  as  well  as  the  surface  wave  contributions.  In  characterizing  these 
waveforms  we  will  prepare  the  ground  for  the  section  on  synthetic  seismograms 
at  far  regional  distances  where  the  goal  will  be  to  realistically  model  the  observed 
wavefield.  If  it  is  possible  to  model  the  shape  of  these  waveforms  even  with 
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relatively  simple  assumptions,  effects  of  source  contributions  to  regional 
seismograms  can  be  quantified. 


2.  Observed  waveforms 

The  observations  used  in  this  study  are  seismograms  from  Nevada  Test  Site 
(NTS)  nuclear  explosions  and  earthquakes  in  the  California/Gulf  of  California 
region.  The  seismograms  were  recorded  at  theLajitas  (LTX)  seismic  station  in  the 
Big  Bend  area  of  southwest  Texas.  One  supplementary  observation  of  an  NTS 
explosion  recorded  at  Albuquerque  is  included.  Station  and  epicenter  locations 
are  shown  in  Fig.l,  which  also  displays  the  great  circle  paths  for  the  waves 
crossing  the  Basin  and  Range  province.  Additional  information  on  each  event  is 
summarized  in  Table  I.  LTX  is  part  of  the  GSETT  experiment  (Golden  &  Herrin 
1989)  and  includes  a  set  of  3-component  velocity  transducers  and  a  set  of  3- 
component  broadband  accelerometers.  The  eigenfrequency  of  the  velocity 
instruments  is  1  Hz,  the  broadband  instruments  are  flat  in  acceleration  response 
between  DC  and  20  Hz.  The  data  were  originally  sampled  at  120  Hz  (SP)  and 
10  Hz  (BB).  For  the  purpose  of  data  reduction  and  consistency  between  short- 
period  and  broadband  seismograms,  the  SP  records  were  subsequently  decimated 
to  a  sampling  rate  of  10  Hz,  as  most  of  the  seismograms  contained  no  energy 
above  5  Hz.  The  effective  bandwidth  for  data  from  both  instruments  extends  to 
3-4  liz  (Fig.2). 

A  total  of  15  events  were  selected  for  the  waveform  analysis  that  follows  (see 
Table  I).  These  events  had  magnitudes  between  4.7  and  5.9,  smaller  events  could 
be  identified,  but  a  signal  to  noise  ratio  of  between  1:2  to  1:3  rendered  them 
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unusable  for  detailed  spectral  analysis  and  modeling.  Six  of  the  recordings  were 
from  nuclear  explosions  at  NTS.  Using  both  earthquake  and  explosion 
recordings  we  hoped  to  minimize  source  bias  in  our  modelling.  In  some 
respects,  however,  waveforms  from  explosions  and  earthquakes  often  showed 
quite  similar  waveform  patterns,  as  will  be  discussed  below. 

Fig.2a  and  b  display  the  vertical  data  for  all  events  recorded  by  SP  (2a)  and  BB  (2b) 
instruments.  The  event  closest  to  the  receiver  is  at  836  km  while  the  most 
distant  event  is  at  2224  km.  Although  sorted  in  ascending  order  with  distance, 
the  records  are  not  displayed  in  a  reduced  time-distance  section,  because  of  the 
uneven  distance  spacing,  with  quite  a  large  number  of  seismograms  in  the  1300- 
1500  km  distance  interval  mainly  from  NTS.  The  first  arrivals  of  the  records 
(either  Pn  (<  1000  km)  or  mantle  P  waves  (>  1000  km)  )  are  aligned  with  a  time 
offset  of  20  sec.  The  record  length  was  chosen  as  600  seconds  in  order  to  include 
all  body  and  surface  wave  contributions  in  the  seismograms. 

As  Koch  &  Stump  (1989)  discussed,  these  v  aveforms  are  characterized  by  an 
emergent  Pn  phase,  strong  Pg  and  Lg  phases  for  distances  less  than  1000  km.  Beyond 
1200  km,  both  P  phases  mentioned  before  can  no  longer  be  identified  and  are 
replaced  by  mantle  P  waves.  This  arrival  along  with  Lg  dominates  the  seismograms. 
Lg  is  further  emphasized  in  this  distance  range  by  the  broadband  recordings,  because 
of  its  lower  frequency  content  relative  to  the  P  waves.  Surface  waves  are  particularly 
strong  for  the  earthquakes,  again  emphasized  in  the  broadband  data  (Fig.2b).  It 
should  also  be  pointed  out  that  there  is  a  lack  of  any  mantle  S  phases.  These  phases 
should  precede  the  Lg  wave,  which  appears  to  be  the  only  identifiable  S  arrival. 
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The  transition  from  crustal  to  mantle  propagation  for  P  waves  is  even  more 
evident  in  fig.3,  where  the  P  wave  segmen  '■'f  a  few  seismograms  is  enlarged.  Pn  is 
very  emergent  for  distances  less  than  1000  km,  and  may  be  seen  as  a  small  arrival 
just  in  front  of  the  mantle  P  wave  for  distances  beyond  1200  km.  In  order  to 
demonstrate  the  frequency  content  of  the  data,  we  calculated  sonograms  from  the 
short-period  data  for  2  events  in  Fig.4.  This  figure  shows  that  the  main  energy  for 
the  regional  seismograms  is  in  the  0.2  to  1.5  Hz  range,  with  a  characteristic  frequency 
decrease  from  P  to  Lg  and  the  surface  waves. 

As  surface  waves  are  often  used  as  discriminants  between  earthquakes  and 
explosions,  we  filtered  the  data  of  Fig.2  with  a  low  pass  filler  at  0.1  Hz  to  enhance 
and  separate  the  surface  wave  contributions.  The  filtered  traces  are  shown  in 
Fig.5  (a,b).  Due  to  the  strong  roll-off  for  the  SP  instruments  below  the 
eigenperiod,  there  is  considerable  high  frequency  energy  remaining  in  the 
seismograms,  while  the  broadband  data  are  characterized  by  the  long-period 
surface  wave  components.  Earthquake  as  well  as  explosion  seismograms  show 
considerable  surface  wave  energy..  Comparing  the  surface  waves  from  both 
sources  in  the  distance  range  between  1336-1468  km,  thjre  are  no  obvious 
differences  that  may  give  rise  to  a  successful  discrimination  between  events. 
Either  the  generation  and  propagation  of  surface  waves  from  explosions  within 
the  Basin  and  Range  is  favorably  enhanced  compared  to  other  geologic  areas,  or 
the  bandwidth  of  the  data  is  not  large  enough  with  respect  to  long  periods  to 
differentiate  between  nuclear  and  tectonic  events.  Most  explosions  showed 
transverse  surface  waves,  that  are  in  phase  with  the  Rayleigh  waves  observed  on 
vertical  and  radial  components  (Koch  &  Stump  1989).  The  Loves  waves  from 
earthquakes  show  higher  phase  velocities,  as  predicted  by  the  dispersion  curves 
derived  by  Priestley  &  Brune  (1978)  from  their  surface  wave  model. 
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The  major  difference  in  the  surface  wave  observations  is  that  the  explosions 
show  considerably  more  higher  frequency  surface  wave  contributions  trailing  the 
fundamental  mode  arrivals  (see  Figs.4,5a).  It  can  be  hypothesized  that  this 
energy  is  trapped  in  shallow  layers  of  the  crust  due  to  the  near-surface  source. 
This  effect  may  be  reinforced  by  the  Western  U.S.  crustal  model  which  is  topped 
by  a  sedin^entary  layer  with  low  seismic  velocities.  Unfortunately  the  signal  to 
noise  ratio  does  not  allow  a  more  detailed  analysis  of  this  pattern. 


3.  Synthetic  seismogram  modelling 

Since  our  goal  is  to  model  the  complete  far-regional  seismograms  with  as  simple 
a  velocity  model  as  possible,  we  used  the  reflectivity  method  (Muller  1985), 
which  allows  us  to  calculate  synthetic  seismograms  including  all  body  and 
surface  waves  for  plane  layered  structures.  We  focused  on  both  amplitude 
information  as  well  as  waveform  shape.  The  velocity  models  were  further 
examined  by  travel  time  calculations  to  facilitate  the  identification  of  particular 
phases  and  to  constrain  the  depth  range  where  arrivals  were  propagating.  For 
the  reflectivity  calculations  throughout  this  study,  we  used  an  explosion  source 
model  with  a  impulsive  time  function,  the  frequency  range  for  the  seismograms 
was  chosen  from  0.05  to  2  Hz,  based  on  the  bandwidth  of  the  data  (Fig.4). 

For  the  Western  U.S  a  number  of  velocity  models  exist  describing  the  structure 
for  the  crust  and  upper  mantle.  These  models  were  derived  primarily  from  long 
period  body  waves  as  well  as  surface  waves,  and,  in  general,  include  the  depth 
range  between  0-500  km.  Some  of  these  models  are  summarized  in  Fig.6,  where 


7 


P  (6a)  and  S  (6b)  models  are  reproduced.  In  order  to  adequately  describe  body  and 
surface  wave  contributions,  we  tried  to  combine  these  models,  as  no  one  by  itself 
was  able  to  satisfactorily  model  our  observed  data.  We  focused  on  the  models  by 
Burdick  and  Helmberger  (1978),  derived  from  P  wave  analysis,  and  Priestley  and 
Brune  (1978),  which  inverted  surface  waves  for  velocity  structure  (P  and  S 
model).  These  two  models  were  previously  combined  and  extended  by  Olsen  et 
al.  (1980)  to  study  wave  propagation  in  the  eastern  Basin  and  Range  for  distances 
between  700  and  900  km.  As  we  started  our  modelling  effort  based  on  these 
models,  a  brief  discussion  of  their  characteristics  is  included. 

The  Burdick  &  Helmberger  (BH)  model  (1978)  consists  of  a  simple  crust 
containing  two  layers  with  P  velocities  of  6.0  and  6.5km/s,  respectively, 
underlain  by  high-velocity  region  below  the  Moho.  This  depth  range,  with  P 
wave  velocities  between  7.95  and  8.05  km/s  has  a  small,  positive  velocity 
gradient.  Below  the  so-called  "mantle  lid"  is  a  low  velocity  layer  with  a  strong 
negative  velocity  jump  to  7.70  km/s,  after  which  the  velocity  increases  slightly 
with  depth.  At  120  -  150  km  a  stronger  velocity  increase  can  be  seen.  Below  this 
depth,  the  velocity  is  again  moderately  increasing.  The  Priestley  and  Brune  (PB) 
model  (1978)  (P  wave  part)  also  has  a  fairly  simple  crustal  structure,  comparable 
to  the  BH  model,  except  for  a  shallow  sedimentary  layer  at  the  surface.  In 
contrast  to  the  BH  model,  the  P  wave  velocities  in  the  PB  model  are  constant 
from  the  Moho  to  a  depth  of  120  km.  Below,  the  velocity  increases  with  a 
moderate  gradient  to  depth  greater  than  300  km.  The  P  wave  part  of  the  PB 
model  is  not  considered  as  an  adequate  model,  ?  ;  the  traveltimes  calailated  from 
this  model  are  too  large  by  more  than  5  sec.  As  Olsen  et  al.  (1980)  pointed  out, 
the  mantle  lid  in  the  BH  model  is  not  well  constrained  by  their  data,  as  stronger 
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mantle  arrivals  for  distances  shorter  than  800  km  are  not  adequately  modeled. 
Therefore,  they  removed  the  mantle  lid  by  a  zone  with  a  small  negative  gradient. 

The  S  wave  part  of  the  PB  model  basically  follows  the  P  velocities  in  the  crust 
with  a  Poisson's  ratio  near  0.25.  Below  the  Moho  a  constant  velocity  layer 
(4.5  km/s)  with  a  thickness  of  about  30  km  follows  with  the  same  P  to  S  ratio.  At 
greater  depth  they  introduced  a  region  with  significantly  lower  velocities  of  4.05 
to  4.10  km/ s.  At  a  depth  of  200  km  a  jump  in  velocity  to  4.4  km/ s  occurs  and  a 
fairly  constant  velocity  gradient  extends  this  model  to  greater  depth. 

Previous  work  has  focused  on  the  analysis  of  P  waves  or  surface  waves  at 
regional  distances.  Few  paper  dealt  with  constraining  the  S  wave  structure  as 
well.  Olsen  et  al.  (1980),  for  example,  used  the  S  wave  model  of  Priestley  and 
Brune  (1978)  in  their  reflectivity  calculations  solely  to  obtain  higher  accuracy  in 
their  P  wave  calculations,  modifying  this  model  in  a  more  intuitive  way.  This 
study  however  tries  to  constrain  the  S  wave  structure  as  well,  in  order  to  gain  a 
more  thorough  insight  into  the  physical  aspects  of  Lg  wave  propagation.  We 
also  tried  to  constrain  some  surface  wave  properties,  although  this  might  be 
hampered  by  the  simple  assumption  of  a  one-dimensional  structure  for  the 
Basin  and  Range.  Surface  waves  are  quite  sensitive  to  the  top  sedimentary 
layers,  which  are  expected  to  vary  greatly  w’lhin  our  distance  range. 

3.1  P  phase 

As  Olsen  et  al.  (1980)  pointed  out,  the  BIl  model  for  P  waves  was  derived  for  NW 
or  SE  propagation  paths  across  the  Basin  &  Range.  As  our  propagation  paths 
were  also  in  the  SE  direction  (Fig.l),  we  calculated  synthetic  seismograms  for 
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both  the  BH/PB  models^  i.e.  using  the  P  wave  model  from  Burdick  &  Helmber- 
ger  (1978)  and  the  S  velocity  model  from  Priestley  &  Brune  (1978),  as  well  as  for 
Olsen  et  al.'s  (1980)  AlO  model.  The  synthetic  seismograms  from  these  calcula¬ 
tions  are  shown  in  Figs.  7  and  8.  We  focus  first  on  the  P  part  of  the  waveforms. 
This  part  of  the  synthetics  for  the  BH  model  shows  very  impulsive  Pn  phases  for 
all  distances  from  its  onset  between  200  and  300  km  to  900  km,  where  Pn  is  the 
first  arrival.  Starting  between  800  and  900  km,  a  secondary  phase  can  be 
identified  coming  out  of  the  upper  mantle..  This  prominent  secondary  phase  is 
not  observed  in  the  Lajitas  data  nor  is  the  Pn  an  impulsive  arrival  (Fig.3). 

The  Olsen  AlO  model  in  contrast  gives  an  emergent  Pn  throughout  the  modeled 
distance  range,  with  Pg  as  the  dominant  phase  out  to  a  distance  of  900  km. 
Beyond  1000  km,  the  mantle  P  phase  becomes  the  strongest  arrival.  This  indica¬ 
tes  that  the  AlO  model  is  better  in  terms  of  the  P  wave  portion  of  our  far-regional 
seismograms  at  Lajitas  for  events  that  propagate  across  the  Basin  and  Range. 

The  Qp  values  used  in  the  reflectivity  calculations  were  selected  for  the  BH 
model  as  200  for  the  upper  crust  and  500-550  for  the  lower  crust  and  upper 
mantle.  For  the  AlO  model  Q  was  taken  according  to  the  figures  given  by  these 
authors.  Therefore  we  selected  200  in  the  upper  crust  and  1000  in  the  lower 
crust/upper  mantle,  except  for  the  depth  range  117-156  km,  where  Olsen  et  al. 
(1980)  introduced  a  Qp  value  of  100.  These  Q  values  and  the  synthetic 
seismograms,  which  were  normalized  to  the  maximum  amplitude  of  each  trace, 
indicate  that  the  observed  waveform  characteristics  are  primarily  dominated  by 
the  velocity  structure  with  minor  effects  from  the  Q  model. 

In  order  to  better  fit  the  observations  we  refined  the  velocity  structure  according 
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to  a  model  that  is  compatible  with  the  AlO  model  right  below  the  Mo  ho  and 
includes  smoother  gradients  to  greater  depths  (Fig.6a).  The  synthetic 
seismograms  resulting  from  this  revised  model  are  presented  in  Fig.  9,  now  only 
for  the  distance  range  from  600-1400  km,  which  is  more  appropriate  for  our 
observation  range.  The  new  model  shows  again  an  emergent  Pn  phase  and  the 
transition  from  crustal  to  upper  mantle  propagation  between  900  and  1000  km. 
This  model  also  gives  a  good  overall  fit  to  the  diminishing  of  Pg  for  distances 
larger  than  1000  km.  Q  values  for  the  crust  were  here  400  while  in  the  upper 
mantle  we  used  a  Qp  between  300  and  350. 

3.2  Constraints  front  S  waves  and  surface  waves 

In  our  effort  to  represent  the  entire  wavefield,  the  models  of  PB  (S  wave  part) 
and  Olsen  et  a.(1980)  do  a  very  poor  job.  In  Fig.  7,  which  used  the  former  S 
velocity  model,  there  is  considerable  energy  arriving  with  a  phase  velocity  of 
4.5km/s,  the  velocity  of  the  high-velocity  layer  underneath  the  Moho.  This 
mantle  phase  is  not  seen  at  all  in  the  observed  data.  Also  the  Lg  phase,  with  a 
velocity  of  about  3.5  km/ s  is  only  a  small  arrival,  although  it  is  the  dominant 
waveform  in  the  observations.  The  AlO  model  with  respect  to  S  velocity  is  only 
slightly  different  from  the  PB  model,  and  therefore  gives  essentially  the  same 
waveform  pattern.  The  Qs  values  used  were  250  throughout  the  crust,  and  450 
in  the  upper  mantle  assuming  values  of  about  half  of  Qp.  For  the  low  Qp  of  the 
asthenospheric  layer,  the  Qs  value  was  decreased  to  75.  The  very  weak  Lg  arrival 
is  therefore  not  primarily  due  to  the  Qs  values  selected,  but  due  to  the  leakage  of 
the  energy  into  the  upper  mantle.  The  crustal  Q  values  are  consistent  with  Q 
values  of  about  200-300  determined  by  Patton  &Taylor  (1984)  from  surface  wave 
analysis  at  about  1  Hz,  as  our  calculations  were  done  for  frequencies  up  to  2  Hz. 


In  order  to  minimize  the  discrepancy  between  the  S  observations  and  the 
reflectivity  synthetic  seismograms,  we  developed  a  new  S  velocity  structure  to 
explain  the  two  basic  features,  seen  for  the  propagation  of  S  waves  across  the 
Basin  &  Range:  (1)  the  dominance  of  the  Lg  phase  over  the  entire  distance  range 
of  the  observations  (800-2000  km);  and  (2)  the  apparent  lack  of  mantle  S  arrivals 
similar  to  the  P  wave  portion  of  the  observed  waveforms.  The  use  of  a 
constrained  S  wave  model  relative  to  the  P  wave  model  (i.e.  for  a  Poisson's  ratio 
of  0.25),  gave  a  result  similar  to  the  other  authors. 

We  perturbed  the  S  velocity  in  the  crust  of  Priestley  &  Brune  (1978)  to  a  constant 
velocity  of  3.5  km/s.  A  high-velocity  lid  below  the  Moho  was  introduced  ,  which 
was  expected  to  act  as  a  strong  reflector  to  avoid  the  leakage  of  energy  out  of  the 
crustal  waveguide.  The  synthetic  seismograms  from  this  model  still  showed 
considerable  mantle  S  energy.  This  S  energy  presumably  is  propagating  in  the 
waveguide  formed  by  the  crust-mantle  boundary  and  the  low  velocity  layer 
beneath  the  high-velocity  lid.  By  systematic  decrease  of  the  thickness  of  the  S  lid 
we  could  decrease  the  mantle  S  wave  contributions  only  to  a  certain  degree.  The 
inclusion  of  a  strong  negative  gradient  at  the  bottom  of  the  lid  as  a  transition  to 
the  low  velocity  layer  resulted  in  a  further  decrease  of  these  mantle  S  arrivals.  In 
addition,  a  low  Qs  value  of  75  was  required  for  the  region  below  the  lid,  with  Qs 
values  inside  the  lid  of  90-100. 

In  order  to  increase  the  duration  of  the  Lg  phase  as  well  as  the  velocity  contrast 
across  the  crust  and  mantle  boundary,  a  laminated  structure  for  S  was  introduced 
within  the  lower  crust.  This  is  similar  to  the  approach  taken  by  Sandmeier  & 
Wenzel  (1986)  to  explain  reverberations  following  crustal  P  wave  arrivals.  The 
synthetic  seismograms  of  this  model  are  shown  in  Fig.9,  where  now  Lg  is  indeed 
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the  anticipated  very  long  lasting  waveform  as  well  as  the  strongest  S  arrival.  The 
weak  arrivals  coming  in  before  Lg  are  the  residual  contributions  from  the  S  wave 
lid  which  might  be  hidden  by  the  P  wave  coda  in  the  observational  data.  These 
arrivals  are  still  large,  if  the  Qs  in  and  below  the  mantle  lid  are  above  100,.  In  this 
case  the  S  waves  leaking  out  of  the  crust  and  the  lid  are  not  satisfactorily 
attenuated. 


4.  Discussion  and  conclusions 

In  order  to  study  the  effects  of  different  sources  on  the  major  regional  phases  we 
have  developed  a  one  dimensional  velocity  model  which  describes  in  both  a 
qualitative  and  quantitative  way  the  propagation  path  effects  for  seismic  waves 
travelling  across  the  Basin  and  Range  from  events  in  Nevada  and  South-Central 
California  to  the  Lajitas  seismic  station  in  the  Big  Bend  of  Texas.  This  model  was 
developed  for  both  P  and  S  broadband  data  in  an  observation  range  larger  than 
used  in  previous  studies  (Olsen  et  al.  1980,  1983).  Previous  models  were  based  on 
long-period  observations,  which  may  be  responsible  for  the  inability  to  resolve 
fine  structures  in  the  crust  and  upper  mantle.  When  comparing  the  BH  and  PB 
models  with  our  new  model,  it  becomes  obvious  that  these  earlier  structures  are 
smoothed  version  of  the  latter.  As  an  example  the  mantle  lid  in  the  BH  model 
was  turned  up-side  down,  i.e.  the  small  positive  velocity  gradient  was  changed 
to  a  small  negative  gradient,  in  order  to  obtain  emergent  Pn  arrivals  throughout 
our  observation  range..  If  the  new  model  is  further  compared  to  the  P  model 
given  by  Priestley  &  Brune  (1978)  it  appears,  that  the  upper  100  km  below  the 
Moho  is  essentially  the  average  of  the  velocities  in  our  model  with  values 
ranging  between  7.7  and  8.0  km/s.  Also,  the  S  velocity  model,  that  these  authors 
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proposed  is  quite  similar  if  our  thin  mantle  S  lid  averaged  over  a  larger  depth 
range.  Furthermore  there  is  a  striking  similarity  of  our  S  model  for  the 
uppermost  mantle  to  the  S  structure  proposed  by  Priestley  &  Brune  (1982)  for  the 
Northern  Volcanic  Plateau.  In  their  comparison  of  different  provinces  within 
the  Basin  and  Range  or  adjacent  areas  they  found  a  thin  mantle  S  lid  with  a 
thickness  of  10  km  and  an  average  S  velocity  of  4.3km/s  correlating  with  our 
results. 

The  main  difference  remaining  between  our  modelled  synthetic  seismograms 
and  the  observations  is  the  amplitude  of  the  Lg  wave.  This  difference  may  result 
from  the  fact  that  the  dominant  source  frequency  in  the  observations  was 
considerably  less  than  2  Hz.  Using  a  source  function  with  lower  frequencies 
would  reduce  the  higher  frequency  P  wave  amplitudes  in  the  synthetic 
seismograms,  but  also  transfer  energy  to  Lg.  This  fact  is  documented  by  the 
seismograms  shown  in  Fig.  10,  where  the  dominant  source  frequency  was 
reduced  from  2  Hz  to  0.5  Hz.  This  is  also  in  accordance  with  the  obseived  Lg  at 
Lajitas,  which  shows,  that  the  main  Lg  contributions  are  near  0.5  Hz  ,  while  the  P 
waves  have  the  largest  amplitudes  in  the  1.0-1 .5  Hz  range  (Fig.4). 

Other  explanations  for  the  strong  Lg  waves  observed  are  related  to  source 
contributions  from  different  secondary  sources.  As  the  synthetic  data  were 
generated  for  a  explosion  source  model,  further  S  wave  contributions  might 
result  from  spall  and  tectonic  stress  release.  Barker  et  al.  (1990)  consider  spall  as 
the  ultimate  source  of  Lg  from  nuclear  explosions.  On  the  other  hand,  our 
calculations  show  that  Lg  can  be  effectively  generated  by  an  isotropic  source  due 
to  conversion  of  the  primary  P  waves.  Further,  the  transfer  of  seismic  energy 
into  Lg  is  more  pronounced  the  lower  the  dominant  source  frequency. 
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Epicentral  distance  to  Lajitas,  TX  (LTX) 


Figure  captions 


Figure  1:  Basin  and  Range  sketch  map  showing  study  area  and  locations  of 
station(s)  and  epicenters  with  corresponding  great  circle  paths. 

Figure  2:  (a)  Short  period  velocity  data  in  the  distance  range  800-2000  km 

(eigenfrequency  of  transducer  - 1  Hz) 

(b)  Broadband  acceleration  data  for  the  same  distance  range 
(bandwidth  0-5  Hz) 

( O  indicates  an  explosion) 

Figure  3:  P  wave  segments  of  selected  events  showing  Pn  and  Pg  for  the  top  2 
traces  and  mantle  P  for  the  lower  traces.  The  weak  precursor  to  the 
mantle  P  waves  is  Pn 

Figure  4:  Sonograms  for  an  earthquake  and  explosion  at  far-regional  distances 
that  illustrate  the  bandwidth  of  different  phases  in  the  seismograms 

Figure  5:  Low-pass  filtered  results  (10  sec)  of  the  data  shown  in  Fig.  2  designed 
to  emphasize  surface  wave  contributions 

Figure  6:  One-dimensional  Basin  &  Range  velocity  models:  (a)  P,  and  (b)  S 
[Burdick  &  Helmberger  1978;  Priestley  &  Brune  1978;  Olsen  et  al.  1980; 
Engen  &  Helmberger  1974] 

Figure  7:  Synthetic  reflectivity  seismograms  at  regional  distances  (200-1400  km) 
for  a  velocity  structure  incorporating  the  P  model  of  Burdick  & 
Helmberger  (1978)  and  the  S  wave  model  of  Priestley  &  Brune  (1978). 


19 


Figure  8: 
Figure  9: 


Figure  10: 


Synthetic  seismograms  for  ve’ocity  models  given  by  Olsen  et  al.  (1980) 

Synthetic  reflectivity  seismograms  at  far-regional  distances  (600- 
1400  km)  for  the  velocity  structure  derived  in  this  study  (see  Fig.  6) 

Same  as  Fig.  9  with  lower  frequency  source  excitation. 
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FIGURE  1 
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Vertical  Component  [BB] 


Vertical  Component  [SP]  Vertical  Component  [BB] 


FIGURE  3 


Frequency  (Hz) 


Sonogram  for  event  at  1336  km 


Sonogram  for  event  at  1468  km 


Vertical  Component  [SP]  /  LP  0.1  Hz 


Vertical  Component  [BB]  /  LP  0.1  Hz 
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P  velocity  (km/s) 


. .  This  study 

Priestley  &  Brune  1978 
.............  Helmberger  &  Engen  1974 

. . .  Olsen  et  at.  1980 

— ' —  Koch  &  Stump  1989 
.  Burdick  &  Helmberger  1978 


FIGURE  6a 


VERTICAL 


figure  7 


VERTICAL 


32 


FIGURE  9 
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Source  Contributions  in  Regional  Waveforms  with  Implications  for 

Explosions  with  Spall 

Karl  Koch,  Brian  Stump 
Southern  Methodist  University 
Department  of  Geological  Sciences 
Dallas,  TX  75275 

1.  Introduction 

A  physical  understanding  of  source  contributions  at  regional  distances  is  an 
important  tool  leading  to  the  identification  and  discrimination  of 
earthquakes  and  explosions  in  the  context  of  a  reduced  threshold  test  ban 
treaty.  With  this  motivation  we  have  analyzed  seismograms  from  nuclear 
explosions  and  earthquakes  at  Lajitas,  TX,  covering  a  distance  range  of  about 
800-2000  km  in  an  attempt  to  explain  observed  waveforms  in  terms  of 
propagation  path  and  source  effects.  Realistic  modelling  of  the  whole 
seismogram  was  carried  out  with  the  reflectivity  method  for  a  specific 
Western  United  States  model  (Koch  &  Stump  1991).  This  work  focused  on 
development  of  both  a  P  and  S  wave  model  that  describes  wave  propagation 
to  far-regional  distances  across  the  Basin  and  Range  for  all  major  regional 
phases.  Based  on  this  work,  we  will  explore  the  role  different  source  models 
play  in  generating  the  regional  phases. 
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The  broadband  data  from  Lajitas  offers  the  opportunity  to  explore  not  only 
the  arrival  times  of  the  dominant  regional  phases  but  their  frequency  content 
as  well.  A  number  of  authors  have  suggested  source  discriminants  based 
upon  the  frequency  content  of  different  regional  phases  (Pomeroy  et  al.,  1982; 
Taylor  et  al.,  1988;  Chael,  1988;  Suteau-Henson  and  Bache,  1988;  Baumgardt 
and  Young,  1990). 

The  initial  deposition  of  seismic  energy  around  a  nuclear  explosion  is 
commonly  represented  as  a  spherically  symmetric  or  isotropic  source.  This 
representation  is  commonly  cast  in  terms  of  a  reduced  displacement  potential 
(rdp)  (Sharpe,  1942;  von  Seggern  and  Blandford,  1972;  and  Mueller  and 
Murphy,  1971).  The  scaling  of  this  rdp  as  a  function  of  yield  is  responsible  for 
the  initial  yield  signature  of  the  seismograms. 

In  near-source  observations,  spall,  the  tensile  failure  of  shallow  layers  caused 
by  interaction  of  the  explosion  pressure  pulse  with  the  free  surface  is  a 
dominant  source  effect.  As  a  secondary  source  of  seismic  waves  it  has  been 
shown  to  have  little  contribution  to  long-period  surface  waves  (Day  et  al. 
1983),  while  studies  of  body  waves  have  been  restricted  to  a  distance  range  of  a 
few  hundred  kilometers  (Taylor  &  Randall  1989,  Barker  et  al.  1990).  In  recent 
papers  it  has  been  suggested  that  this  secondary  source  may  be  responsible  for 
the  success  of  some  spectral  discriminants  (Taylor  &  Randall  1989)  at  regional 
distances.  This  secondary  source  has  also  been  suggested  as  a  strong 
contributor  to  the  Lg  phase  in  some  geologies  (Barker  et  al.  1990). 

In  order  to  quantify  source  effects  at  far-regional  distances,  i.e  distances  of 
more  than  1000  km,  we  have  calculated  complete  seismograms  from 
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explosion  and  spall  sources  for  a  variety  of  different  source  characteristics. 
First  we  want  to  contrast  spall  and  explosion  wavefields  in  terms  of  source 
depth,  which  may  be  different  for  a  normally  contained  explosion  and  the 
near-surface  spall  phenomena.  In  this  parametric  study  we  use  the  same 
source  time  functions  in  order  to  separate  depth  effects  from  effects  that  may 
result  from  spectral  source  differences. 

Near-source  observations  from  two  nuclear  events  at  Rainier  Mesa,  NTS, 
will  be  incorporated  in  the  modelling  of  far-regional  seismograms  for  a  range 
of  source  parameters.  Both  free-field  and  surface  observations  "  /ere  available, 
which  we  use  to  constrain  the  source  ’'irameters  of  both  explosion  and  spall 
components  in  terms  of  a  von  Seggern-Blandford  explosion  source  model 
(von  Seggern  &  Blandford  1972)  and  a  spall  source  model  according  to  Stump 
(1985). 

2.  Spectral  analysis  of  observed  waveforms 

In  order  to  quantify  spectral  effects  that  may  be  identified  in  regional 
seismograms  we  analyzed  the  data  from  the  Lajitas,  TX  seismic  station,  which 
are  shown  in  Fig.l.  These  data  were  recorded  by  a  broadband  acceleration 
instrument  with  a  bandwidth  from  DC  to  20  Hz  but  noise  levels  and  anti-alias 
filters  limit  the  bandwidth  of  the  data  to  3-4  Hz.  These  data  are  a  subset  of 
those  used  by  Koch  &  Stump  (199i)  to  constrain  propagation  path  effects  at 
far-regional  distances  across  the  Basin  &  Range.  They  are  distinguished  by 
signal  to  noise  ratios  of  at  least  5.  Both  explosion  and  earthquake  data  were 
used  for  the  spectral  analysis  in  order  to  discuss  spectral  differences.  Three 
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explosion  recordings,  marked  by  an  asterisk  in  Fig.  1,  were  compared  to  five 
earthquake  recordings. 

To  characterize  the  temporal  development  of  the  spectral  content  in 
explosion  and  earthquake  seismograms  we  calculated  sonograms  from  the 
raw  data.  A  time  window  of  2.56  sec  was  used  for  each  individual  spectral 
estimate,  and  for  smoothing  purposes,  this  window  was  moved  in  time  by 
0.64  sec  for  subsequent  spectra,  resulting  in  a  fourfold  oversampling  and 
smoothing  of  each  sonogram.  The  sonograms  were  plotted  in  3-dimensional 
perspective  view  and  as  contour  plots.  Representative  figures  for  two  events 
used  in  our  analysis  are  shown  in  Figs.2  and  3,  where  the  first  is  from  a  NTS 
explosion  at  a  distance  of  1468  km  ("Kearsarge"  -  referred  to  in  the  following 
as  trace  no.6)  and  the  second  event  is  an  earthquake  (trace  no.7)  in  Southern 
California  at  a  similar  range  (1504  km).  Due  to  the  proximity  in  epicentral 
distance,  wave  propagation  effects  are  assumed  to  be  similar  allowing  a  direct 
comparison  for  these  events. 

The  sonogram  for  event  #6  (Fig.2)  shows  distinct  frequency  signatures  for 
each  regional  phase  including  mantle  P,  Lg  and  surface  waves.  The  early  part 
of  the  sonogram  is  dominated  by  the  mantle  P  wave  with  energy  in  the 
frequency  band  0.3  to  1.5  Hz.  The  Lg  arrival  occurs  at  about  220-230  sec  with  a 
*  shift  of  the  energy  to  lower  frequencies,  and  a  strong  spectral  maximum 

around  0.5  Hz.  At  times  greater  than  300  sec  surface  waves  arrive  with 
further  frequency  shift  to  0.1-0.2  Hz.  As  the  seismograms  in  Fig.l  suggest,  the 
sonogram  of  the  earthquake  record  (event  #7)  is  strongly  dominated  by  the 
surface  wave  arrivals.  This  arrival  begins  after  300  sec  and  the  frequency 
content  of  these  wave  is  similar  to  those  from  the  explosion.  The  Lg  phase. 


37 


starting  at  250  sec,  again  shows  a  dominant  frequency  range  around  0.5  sec. 
The  P  waves,  due  to  the  dominance  of  the  surface  waves,  are  very  small  in 
the  sonogram  although  they  contain  the  highest  frequencies,  up  to  1  Hz.  This 
high  frequency  phase  characterization  is  further  supported  by  sonograms  of 
short-period  data  recorded  from  these  events  which  suppress  the  surface 
waves  below  0.3-0.5  Hz. 

Taylor  et  al.  (1988)  and  Chael  (1988)  have  suggested  that  the  spectral  shape  for 
earthquake  and  explosion  sources  is  different.  This  spectral  discriminant  was 
tested  by  calculating  the  P  wave  spectra  of  the  explosions  (Fig.4a)  and 
earthquakes  (Fig.4b)  from  the  broadband  seismograms.  A  time  window  of  80  - 
150  sec  was  used  in  the  spectral  calculation  incorporating  alt  P  wave 
contributions  -  Pn,  Pg,  mantle  P.  No  attempt  was  made  to  separate  these 
phases  due  to  the  following  reasons:  (1)  the  data  base  currently  available  is  too 
small  for  a  useful  phase  separation,  i.e.  only  two  seismograms  show 
reasonable  signal  levels  in  Pn  and  Pg,  while  all  other  events  show  only 
mantle  P  waves;  (2)  as  the  focus  of  this  work  is  on  source  effects,  source 
contributions  should  be  contained  in  all  P  phases  and  hence  should  also  be 
observed  in  this  integral  measure.  As  we  used  broadband  acceleration  data, 
the  spectra  were  corrected  for  instrument  effects  by  double  integration  in  time 
transforming  them  into  displacement  spectra.  In  order  to  separate  the 
individual  event  spectra  in  the  figures,  they  were  scaled  by  multiple  factors  of 
10.  The  spectra  were  not  corrected  for  epicentral  distance  or  anelastic 
attenuation,  nor  were  they  normalized  with  respect  to  source  strength. 

While  the  earthquakes  (Fig.4a)  show  a  clear  spectral  plateau  for  low 
frequencies,  the  explosions  (Fig.4b)  are  characterized  by  a  decrease  of  spectral 
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amplitudes  for  low  frequencies,  in  agreement  with  the  explosion  source 
models  of  Mueller  &  Murphy  (1971)  or  von  Seggern  &  Blandford  (1972).  To 
further  emphasize  the  spectral  differences  we  used  the  closely  spaced  events 
in  the  1300-1500  km  distance  range  for  spectral  comparison  (Fig.5).  The 
decrease  in  long-period  spectral  amplitude  for  the  explosions  in  contrast  to 
the  earthquakes  is  obvious.  The  corner  frequencies  of  the  explosions 
apparently  are  shifted  to  higher  frequencies  and  the  high-frequency  roll-off  is 
larger.  As  Taylor  et  al.  (1988)  argued,  this  difference  in  spectral  amplitudes 
might  be  used  as  a  discriminant,  and  they  proposed  the  spectral  ratio  in  the  1- 
2  Hz  and  6-8  Hz  bands.  The  spectra  for  our  data  indicate  such  a  measure  may 
be  useful  at  far-regional  distances  for  the  frequency  bands  at  0.05-0.2  Hz  (5-20 
sec)  and  0.3-1  Hz,  where  the  explosions  show  the  largest  difference  compared 
to  the  plateau  for  earthquakes. 


3.  Source  depth  effects  for  explosion  and  spall  sources 

Near-source  ground  motions  (ranges  to  several  kilometers)  from  nuclear 
explosions  are  dominated  by  the  explosive  source  itself  and  the  effects  of 
spallation  of  near-surface  layers.  In  order  to  quantify  these  spall  effects.  Day  et 
al.  (1983)  developed  a  model  to  describe  the  tensile  failure  of  near  surface 
layers.  Their  equivalent  source  time  function  is  composed  of  the  initial 
impulsive  failure  of  the  near-surface  layers,  the  relaxation  during  frep  fall  of 
the  spalled  layers,  and  the  slapdown  pulse  corresponding  to  the  rejoin  of  the 
spalled  mass.  The  driving  mechanism  is  usually  either  described  by  a  single 
vertical  force  (f3)  or,  often  applied  in  moment  tensor  inversions,  by  the 
equivalent  force  couple  model  (mn,  m22, 1^33)  (Day  &  McLaughlin  1991).  In 
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our  modeling  of  far-regional  seismograms  we  have  adopted  the  force  couple 
model  for  the  spall  source.  The  source  time  functions,  Mueller-Murphy  (1971) 
or  von  Seggern-Blandford  (1972)  for  the  explosion  source,  and  Day  et  al.(1973) 
or  Stump  (1985)  for  the  spall  source,  were  chosen  as  a  simple  step  function  in 
order  to  focus  on  source  depth  effects  alone.  We  will  consider  specific  source 
time  functions,  along  with  particular  source  depths  for  primary  explosion 
source  and  secondary  spall  source  in  a  following  section. 

In  order  to  quantify  the  depth  dependence,  we  calculated  synthetic  reflectivity 
seismograms  for  a  variety  of  source  depths  from  0-2  km.  The  depths  for 
which  synthetic  seismogram  calculations  were  performed  included  0  m,  100 
m,  300  m,  500  m,  750  m,  1  km  and  2  km  (Fig.6).  Shown  in  this  figure  is  the 
near-source  velocity  model  consisting  of  a  low  velocity  layer  (2.5  km  thick) 
with  P  and  S  wave  velocities  of  3.5  km/ s  and  2.09  km/s,  followed  by  a  typical 
crustal  velocity  structure  with  the  corresponding  velocities  of  6  km/s  and  3.55 
km/ s.  All  sources  were  placed  in  the  top  low  velocity  layer  typical  of  a  Basin 
&  Range  model  and  the  NTS  area.  Besides  the  source  model,  the  primary 
source  effects  are  expected  to  result  from  the  interaction  of  the  wavefield  with 
the  free  surface  and  the  interface  between  the  low  velocity  sedimentary  layer 
and  crustal  layer.  The  explosion  source  was  modelled  by  an  isotropic  moment 
tensor.  The  spall  source,  physically  describing  the  tensile  failure  of  the 
shallow  layers,  i.e.  the  opening  of  a  horizontal  crack,  is  theoretically 
represented  by  an  isotropic  moment  tensor  with  an  additional  m3  3 
component.  The  ratio  between  the  vertical  and  the  horizontal  force  couples  is 
of  the  order  of  2-3.  As  our  interest  is  in  the  difference  between  the  explosion 
and  spall  event,  we  removed  the  isotropic  component  from  the  spall  source 
model,  using  only  the  m33  component.  We  therefore  assume  that  linear 
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superposition  of  mas  with  the  isotropic  part  would  be  appropriate,  but  that  it 
provides  no  further  information  for  discriminating  between  the  spall  and 
explosion  source.  Differences  in  these  calculations  for  an  explosion  and  our 
spall  source  may  thus  represent  an  upper  bound  for  differences  observed  in 
regional  seismograms. 

From  the  synthetic  data  in  the  distance  range  200-1400  km,  we  determined  the 
maximum  amplitudes  for  the  major  regional  phases.  In  order  to  quantify 
depth  and  range  effects  in  the  synthetic  seismogram  sections  for  the  various 
wave  types,  amplitude  distance  curves  were  derived  from  the  data.  Fig.7a 
contains  amplitude  distance  curves  for  the  vertical  motion  from  explosion 
sources.  The  amplitude  distance  curves  for  LR  (surface  wave)  show  almost 
identical  values  for  sources  in  the  depth  range  between  0  m  and  2  km,  while 
the  variability  of  the  Lg  and  Pg  amplitudes  are  a  factor  of  2-3,  with  amplitudes 
increasing  with  source  depth.  The  2  km  deep  source  shows  significantly 
higher  amplitudes  for  the  body  waves  which  might  be  expected  from  the 
proximity  to  the  underlying  high  velocity  layers  (see  Fig.6).  Pn  and  mantle  P 
wave  (Pm)  are  plotted  as  one  curve,  and  they  indicate  the  transition  between 
these  two  waves  as  first  arrival.  In  addition,  these  waves  show  a  substantial 
depth  dependency  as  a  result  of  increased  coupling  for  the  deeper  sources.  The 
Pn/Pm  curves  clearly  indicate  that  the  transition  of  Pn  to  Pm  as  the 
prominent  first  arrival  occurs  between  900  and  1000  km,  independent  of 
source  depth. 

Amplitude  distance  curves  for  vertical  component  seismograms  from  the 
spall  source  are  given  in  Fig.7b.  The  change  in  synthetic  amplitudes  with 
source  depth  is  much  greater  than  those  observed  from  the  explosion  source. 


This  fact  is  emphasized  by  the  amplitude  distance  curves,  which  show 
strongly  increasing  maximum  amplitudes  as  the  source  is  moved  away  from 
the  free  surface.  The  amplitude  variations  with  respect  to  the  spall  source 
depth  is  of  the  order  of  10  compared  to  a  factor  of  2-3  for  the  explosion 
sources.  Due  to  wavenumber  filtering  in  the  synthetic  seismograms,  the 
curves  for  LR  are  not  well  constrained. 

To  further  emphasize  source  differences,  the  vertical  component  amplitude 
distance  curves  for  Pg  as  well  as  Lg  for  explosion  and  spall  sources  are 
superimposed  in  Fig.8a,b.  Scaling  of  the  amplitudes  between  the  two  source 
types  is  relative.  Contrasting  Pg  amplitudes  in  Fig.Sa,  the  differences  for  the 
explosion  source  are  small  except  for  the  very  deep  source  while  the  spall 
source  differences  are  as  large  as  a  factor  of  10.  The  same  comparison  is  found 
for  the  Lg  amplitudes  in  Fig. 8b.  These  results  indicate  that  body  wave 
amplitudes  for  a  deep  spall  source  could  be  quite  large  in  far-regional  data 
while  a  shallow  spall  source  may  not  be  identifiable  in  these  waveforms.  As 
Stump  (1985)  showed,  the  spall  source  has  its  largest  depth  extend  at  ground 
zero  and  may  be  on  the  order  of  one  half  the  burial  depth.  As  the  seismic 
moment,  and  hence  source  strength,  are  determined  by  the  spalled  mass  and 
its  escape  velocity,  the  center  of  the  spall  zone  may  make  a  large  contribution 
to  the  generation  of  significant  seismic  energy  observable  at  regional 
distances. 

4.  Misty  Echo  /  Mineral  Quarry:  A  case  example  for  near-source  constraints 
from  nuclear  explosion 
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Modeling  of  regional  data  alone  introduces  ambiguities  in  the  relative 
importance  of  the  different  source  contributions  such  as  the  isotropic 
explosion  and  the  secondary  spall  source.  In  an  attempt  to  further  quantify 
the  source  function  for  the  regional  seismograms,  near-source  observations 
can  be  used  to  constrain  source  parameters  of  the  explosion  and  spall  sources. 
Combined  free-field  (shot  level)  and  surface  data  further  provides  separation 
of  the  different  source  models.  Patton  (1990),  for  example,  has  used  surface 
data  to  constrain  spall  source  parameters  for  Pahute  Mesa  explosions.  A 
similar  approach  will  be  used  here  for  the  two  nuclear  explosions  Misty  Echo 
and  Mineral  Quarry  at  Rainier  Mesa.  Data  from  these  events  provide  us  with 
the  opportunity  for  quantifying  the  isotropic  source  with  free-field  data  and 
the  spall  source  with  free  surface  observations  from  within  the  spall  zone. 

The  instrumentation  plan  for  Misty  Echo  and  Mineral  Quarry  are  shown  in 
Fig.9  a  and  b.  For  Misty  Echo  there  were  8  acceleration  gauges  in  the  Rainier 
Mesa  tunnel  complex  providing  free-field  data,  in  addition  to  18  surface 
installations,  which  were  used  to  map  the  spall  zone.  For  Mineral  Quarry,  15 
instruments  were  emplaced  in  the  tunnels,  and  more  than  30  surface  stations 
were  operated. 

The  explosion  source  parameters  were  estimated  by  Min  &  Stump  (1991) 
using  a  simultaneous  inversion  of  near-source  data  for  source  parameters 
and  Q.  The  source  was  either  constrained  to  a  Brune's  or  a  von  Seggern  & 
Blandford  model.  While  the  former  model  was  developed  for  earthquakes,  it 
represents  the  limiting  case  of  a  von  Seggern  and  Blandford  model  with  no 
overshoot.  The  von  Seggern-Blandford  model  is  usually  parameterized  by 
two  values  B  and  k,  where  k=2ji-fc  is  the  angular  corner  frequency  of  the 
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source  spectrum,  and  B=(2  A+1)  is  the  overshoot  ratio.  Min  &  Stump  (1991) 
estimate  a  corner  frequency  of  about  1  Hz  for  Misty  Echo  while  the  corner 
frequency  for  Mineral  Quarry  is  1.5  Hz,  consistent  for  both  a  Brune  and  von 
Seggern  &  Blandford  model.  The  overshoot  value,  B,  was  not  as  well 
constrained,  as  it  trades  off  with  the  long-period  level  in  their  inversions. 
Reasonable  values  for  B  between  2  and  5  were  found.  In  our  parametric  study 
of  explosion  and  spall  synthetics  in  the  following  section,  we  use  a  range  of 
1-2  Hz  for  fc  and  1-5  for  B. 


Free  surface  acceleration  records  from  directly  above  the  explosion  constrain 
the  spall  source  parameters.  These  data  were  converted  to  velocity  in  order  to 
determine  an  escape  velocity  estimate  for  the  spalled  mass.  Fig.lO  shows  a 
typical  example  of  the  surface  data,  where  the  station  was  within  the  spall 
zone.  It  exhibits  the  characteristic  -Ig  dwell  during  the  free  fall  of  the  spalled 
mass.  Rise  time,  i.e.  the  time  for  the  spalled  mass  to  detach,  and  dwell  time, 
i.e.  the  time  of  free  fall,  were  also  obtained  from  the  acceleration  data.  The 
dwell  time  was  compared  to  the  escape  velocity  for  consistency: 


Td  = 


2Vo 

g 


where  Tp)  is  the  dwell  time,  Vq  is  the  escape  velocity,  and  g  is  gravitational 
force.: 


Fig.  11  compares  the  vertical  velocities  from  both  free-field  and  free  surface 
data.  In  general,  the  velocities  from  the  free-field  data  for  Misty  Echo  are  a 
factor  2  larger  than  the  values  observed  for  Mineral  Quarry.  In  contrast,  the 
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surface  velocities  inside  as  well  as  outside  the  spall  zone  are  almost  identical 
for  both  events.  The  spall  radius,  however,  determined  from  the  data,  shows 
again  a  factor  2  difference,  where  we  estimate  a  spall  radius  of  1000  m  for 
Misty  Echo  and  500  m  for  Mineral  Quarry.  The  peak  velocities  shown  in  Fig. 
12  indicate  a  peak  escape  velocity  of  0.2-0.3  g-s  for  Misty  Echo  and  values  of 
0.15-0.2  g-s  for  Mineral  Quarry.  In  order  to  bound  spall  momentum  and 
determine  if  this  secondary  source  can  significantly  contribute  to  the  regional 
waveforms  a  range  of  momentum  estimates  for  the  process  were  made.  The 
upper  bound  of  these  values  assumes  that  the  entire  spall  zone  moves  with 
the  escape  velocity  determined  at  the  center  of  the  spall  zone. 

The  dwell  times  are  shown  as  a  function  of  free  surface  range  for  Misty  Echo 
in  Fig.l3.  From  these  figures  we  estimate  a  maximum  dwell  time  between 
0.75  and  0.5  s  for  Misty  Echo.  The  dwell  times  for  Mineral  Quarry,  0.375  s,  are 
about  a  half  as  long  as  those  for  Misty  Echo.  From  these  data  our  estimates  for 
the  average  escape  velocity  are  3.7  m/s  for  Misty  Echo,  and  1.9  m/s  for 
Mineral  Quarry.  Both  of  these  estimates  are  greater  than  the  actual  escape 
velocities  observed  and  indicate  that  the  spall  process  may  not  be  totally 
explained  by  the  simple  ballistic  model  linearly  linking  escape  velocity  and 
dwell  time. 

The  values  independently  derived  are  now  used  to  calculate  estimates  of  the 
spalled  mass  and  the  total  momentum.  They  are  contrasted  against  values, 
that  follow  from  scaling  relations  given  by  Patton  (1990).  His  scaling  relations 
give  estimates  for  spall  extent  as  well  as  depth.  The  maximum  range  of  spall 
is  given  by: 

Tmax  =  475(±60)  W  (±0-03) 
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and  the  depth: 

dmax  =  86(±12)  W  0-25  (±0-03) 

Following  the  suggestion  of  Patton  (1990),  we  estimate  the  spall  mass  by  two 
discs  with  lateral  extent  of  r^ax  and  rmax/2/  respectively,  and  a  depth  of  each 
disc  of  one-half  the  depth  extent  determined  for  the  entire  spall  zone,  hence 

M  =  itpir2d5  +7tp2r^d^  + 

and  we  compare  it  with  his  empirical  relation  for  mass: 

M  =  7.3  X  10^0  w  0.77  kg 
for  shots  below  the  water  table. 

The  results  for  the  two  nuclear  explosions  under  investigation  are 
summarized  in  Table  1.  The  scaling  relations  are  based  on  Lg  yields  reported 
oy  Patton.  The  YLg  for  Mineral  Quarry  is  taken  to  be  6.1  ktons  and  for  Misty 
Echo,  16.9  ktons.  In  the  case  of  the  bigger  of  the  two  events.  Misty  Echo,  the 
total  mass  as  predicted  by  the  scaling  relations  is  within  7%  of  that 
determined  by  the  spall  zone  data.  The  Mineral  Quarry  spall  zone  data 
indicates  a  unusually  small  mass  in  comparison  to  the  scaling  relation  where 
the  difference  is  70%.  This  difference  is  despite  using  a  large  estimate  for  spall 
depth  in  the  data  estimates  which  was  set  to  1/2  the  depth  of  burial  for  the 
explosion. 

5.  Misty  Echo;  Regional  observations  and  synthetic  seismograms 
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As  the  work  presented  in  previous  sections  indicated,  spall  as  a  secondary 
source  may  be  important  and,  in  some  way,  contribute  to  regional 
seismograms.  Barker  et  al.  (1990)  used  a  numerical  simulation  approach  for 
explosion  and  spall  sources  and  found  that  spall  can,  in  some  instances,  be  a 
governing  source  for  the  generation  of  Lg  waves  observed  at  regional 
distances.  While  spall  was  an  efficient  generator  of  Lg  waves  for  a  high- 
velocity  near-surface  structure,  similar  to  what  is  found  for  the  Eastern 
Kazakhstan  test  site,  explosion  sources  were  poor  generators  of  Lg  waves, 
although  capable  of  producing  the  P  wave  portion  of  regional  and  teleseismic 
seismograms.  However,  for  a  velocity  structure  with  a  low  velocity  surface 
layer,  their  model  calculations  show  that  Lg  generation  at  regional  distances 
can  be  attributed  to  either  the  explosion  or  spall  source.  Our  goal  is  to  use  the 
near  source  information  gained  for  the  explosion  and  spall  source  models  to 
calculate  synthetic  reflectivity  seismograms  at  far-regional  distances  and 
investigate  their  significance  on  regional  seismograms. 

Far-regional  observations  from  the  Misty  Echo  explosion  at  Lajitas  seismic 
station  are  shown  in  Fig.l4,  with  a  station-event  distance  of  1455  km.  All 
three  component,  short-period  seismograms  are  shown.  As  NTS  is  at  a 
backazimuth  of  aoout  290  degrees  from  this  station,  the  EW  component  is 
primarily  radial,  while  NS  is  nearly  transverse.  The  data  indicate  a  well 
developed  mantle  P  phase  as  well  as  Lg,  which  is  the  biggest  arrival.  What 
might  be  anticipated  from  the  discussion  in  section  2  and  by  the  seismograms, 
is  a  shift  to  lower  frequency  from  P  waves  to  Lg. 

The  source  functions  used  in  the  synthetic  seismograms,  as  derived  from 
near  source  data,  are  given  in  Fig.l5.  The  explosion  is  represented  in  terms  of 


a  von  Seggern-Blandford  explosion  model  (top)  and  the  spall  source  (bottom) 
follows  the  representation  of  Stump  (1985).  For  the  explosion  source  model, 
the  far-field  displacement,  or  equivalent  near-field  velocity,  is  given  for  a 
number  of  k  values,  corresponding  to  corner  frequencies  of  1,  1.5  and  2  Hz. 
The  source  time  functions  are  scaled  to  a  normalized  static  reduced 
displacement  potential.  For  the  spall  source  we  have  plotted  the  ground 
acceleration  in  g's  while  dwell  and  rise  times  vary  between  0.375  and  0.75  s. 
With  these  source  parameters  we  calculated  regional  synthetics  for  the 
distance  range  from  600  to  1400  km.  For  emplacement  depth  we  used  500  m 
for  the  explosion,  close  to  the  actual  source  depth,  while  we  used  200  km  for 
the  spall  source,  which  we  have  estimated  from  the  near  source  data. 

Fig.l6  shows  reflectivity  seismograms  for  a  explosion  source  with  no 
overshoot  (B=l)  and  a  dominant  frequency  around  1  Hz  (k=6),  which 
correspond  to  the  results  obtained  from  Koch  &  Stump  (1991)  for  their  study 
of  wave  propagation  effects  across  the  Basin  and  Range.  As  no  specific  source 
information  was  incorporated,  this  figure  illustrates  the  partitioning  of 
seismic  energy  into  different  regional  phases  that  might  be  enhanced  or 
blurred  by  source  depth,  source  mechanism,  and  source  time  functions.  A 
further  result  provided  by  these  synthetic  data  is  the  higher  frequency  content 
for  Pn,  Pg  and  mantle  P  waves  when  compared  to  the  Lg  and  surface  wave 
parts  of  the  seismograms. 

Fig.l7  shows  a  similar  synthetic  seismogram  section,  where  the  overshoot 
ratio  was  increased  to  B=5,  thus  strongly  enhancing  the  energy  around  1  Hz. 
As  the  synthetics  show  this  results  in  a  further  enrichment  of  Pg  and  mantle 
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P  energy,  while  the  Lg  as  well  as  the  surface  waves  remain  almost  the  same, 
due  to  their  longer  period  nature. 

For  the  spall  source  we  show  the  results  for  the  distance  of  1400  km  with  a 
parametric  change  of  the  dwell  and  rise  time  in  the  source  time  functions. 
Our  estimates  for  these  two  parameters  were  between  0.375  and  0.5  sec  for  the 
dwell  time  and  0.5-0.75  sec  for  the  rise  time.  We  also  used  upper  limit  values 
twice  as  large  to  bound  the  effects  on  the  regional  waveforms.  Fig.18 
illustrates  the  strong  dependency  of  the  seismograms  on  these  source 
parameters.  Although  only  a  qualitative  examination  is  made,  the  P  wave 
appears  to  be  a  more  stable  contribution  to  regional  waveforms  than  either 
the  Lg  wave  or  surface  waves.  This  may  be  related  to  the  low  frequency 
content  of  these  later  phases.  The  spall  source  function  is  peaked  and  moving 
this  peak  to  higher  frequencies  strongly  affects  the  later  phases.  Barker  et  al. 
(1990)  argue  that  the  peak  of  the  spectrum  for  the  spall  source  corresponds  to 
the  dwell  time,  In  this  light,  spall  might  not  be  efficiently  contributing  to 
lower  frequency  Lg  or  surface  waves  except  in  the  case  of  large  dwell  times 


6.  Conclusions  and  Summary 


In  order  to  quantify  the  effects  of  spall  on  regional  seismograms,  which  in  our 
case  are  far-regional  distances  from  several  hundred  to  about  2000  km,  we 
investigated  data  from  the  Lajitas  seismic  station.  We  conclude  from  our 
spectral  analysis,  that  regional  P  waves  across  the  Basin  and  Range  show 
significantly  higher  frequencies  than  Lg  waves,  although  the  latter  appear  as 
strong  arrivals  especially  in  short  period  data. 

As  recent  work  (Barker  et  al.  1990)  identified  spall  as  a  possible  strong 
generator  of  Lg  waves  at  regional  phases,  this  study  has  focused  on  the 
explosion  and  spall  parameters  that  might  separate  their  contributions  to 
regional  seismograms  and  quantify  their  significance.  We  first  tried  to  tie 
down  the  significance  of  source  depth  effects  from  both  source  types.  As  was 
demonstrated,  there  is  a  strong  depth  dependency  for  excitation  of  all  regional 
phases  due  to  the  spall  source,  while  the  differences  for  the  explosion  are  less 
distinct. 

As  spall  is  a  major  contribution  to  near  source  observations,  we  used 
observations  of  the  nuclear  explosions  Misty  Echo  and  Mineral  Quarry  to 
constrain  explosion  and  spall  source  parameters  for  small  to  medium  sized 
events.  Using  these  constraints  we  observe  that  spall  may  not  necessarily  be  a 
significant  contribution  to  all  regional  seismic  waves.  We  attribute  ^his  result 
to  the  low  frequency  character  of  observed  Lg  waves,  which  may  not  be 
generated  by  spall  sources  with  short  duration  and  vanishing  momentum  at 
longer  periods. 
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Figure  Captions 


Figure  1:  Far-regional  seismograms  (broadband  acceleration  data)  for  the 
distance  range  800-1500  km;  records  from  explosions  are  marked  by 
the  O  symbol. 

Figure  2:  Soi\ogram  for  an  NTS  explosion  at  a  distance  of  1468  km;  (top) 
perspective  view,  (bottom)  contour  plot.  Time  and  frequency  axes 
have  the  same  scale  in  both  plots. 

Figure  3:  Sonogram  for  an  earthquake  in  Southern  California  at  a  distance  of 
1504  km;  same  format  as  Fig.2. 

Figure  4:  P  wave  spectra  from  (a)  earthquake  and  (b)  explosion  seismograms. 

A  window  of  80-150  sec  has  been  used  in  order  to  include  the 
complete  P  wave.  Spectra  are  offset  by  a  factor  of  10  for  display 
purposes. 

Figure  5:  Comparison  of  P  wave  spectra  for  explosions  and  earthquakes  at  a 
range  of  1400-1500  km. 

Figure  6:  Source  and  emplacement  model  used  for  investigating  source  depth 
effects  between  0  and  2  km. 

Figure  7:  Amplitude  distance  curves  for  (a)  explosion  and  (b)  spall  sources  in 
the  distance  range  200-1400  km. 
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Figure  8:  Comparison  of  Pg  and  Lg  amplitude  distance  curves  for  different 
source  models. 

Figure  9:  Map  of  instrument  setup  for  near  source  survey  for  (a)  Misty  Echo 
and  (b)  Mineral  Quarry. 

Figure  10:  Near  surface  acceleration  recording  and  corresponding  velocity 
seismogram  calculated  by  integration. 

Figure  11:  Peak  vertical  velocities  versus  slant  range  for  Misty  Echo  and 
Mineral  Quarry  using  both  free-field  and  surface  data. 

Figure  12:  Peak  vertical  velocities  versus  free  surface  range  for  free  surface 
data.  The  plots  were  used  to  constrain  escape  velocity  and  spall 
dimension. 

Figure  13:  Dwell  times  for  spall  region  data  from  Misty  Echo  and  Mineral 
Quarry 

Figure  14:  Three  component  far-regional  seismogram  at  Lajitas  (LTX)  from 
the  NTS  explosion  Misty  Echo;  epicentral  distance  1455  km. 

Figure  15:  Source  time  functions  constrained  by  near-source  observations;  a 
von  Seggern-Blandford  model  was  adopted  for  the  explosion  source 
with  B=5  and  k=6,  9  and  12;  for  the  spall  source  the  model  of  Stump 
(1985)  was  used  with  dwell  times  ranging  from  0.5-0.75  sec,  and  rise 
times  from  0.375-0.75  sec. 
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Figure  16:  Far-regional  synthetic  seismograms  for  the  velocity  model 
developed  by  Koch  &  Stump  (1991)  from  an  explosion  source 
without  overshoot. 

Figure  17:  Far-regional  synthetic  seismograms  for  a  von  Seggern-Blandford 
explosion  source  with  k=6  and  B=5. 

Figure  18:  Far  regional  seismograms  for  spall  sources  with  different  source 
parameters  (dwell  time,  rise  time)  at  a  distance  of  1400  km. 
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Total  Momentum:  2.9-5.5  x  10^2  Nt- 


Vertical  Component  [BB] 


57 


FIGURE  1 


Event  #  6  [Component=Z-BB] 


100 


Time  (sec) 
FIGURE  2 
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Event  #  7  [Component=Z-BB] 


Time  (sec) 
FIGURE  3 
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1896  km 


Frequency  (Hz) 

FIGURE  4a 


1468  km 


Frequency  (Hz) 


Synthetic  Seismograms  *  Explosion  Source 
Vertical  Component  (0.05  -  2  Hz) 
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Synthetic  Seismograms  *  Spall  Source 
Vertical  Component  (0.05  -  2  Hz) 
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FIGURE  8a 


Comparison  of  Amplitudes  for  Explosion  and  Spall  Source 
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FIGURE  8b 
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FIGURE  11 
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MISTY  ECHO 


Time  (sec) 

FIGURE  15a 


MISTY  ECHO 
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FIGURE  15b 


200  400  s 

FIGURE  18 
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TELESEISMIC  CODA  ANALYSIS 


WILLIAM  L.  SOROKA 
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PAPER  I 


SUMMARY 

Spectra  from  teleseismic  coda  waves  are  used  to 
estimate  earth  and  source  properties  through  Generalized 
Linear  Inversion  (GLI)  techniques.  High  frequency 
teleseismic  (66-98°)  P-wave  codas  recorded  by  the  Regional 
Seismic  Test  Network  (RSTN)  are  analyzed.  The  spectra  of 
the  teleseismic  coda  were  found  to  be  described  by  the 
following  two  models:  (A/t^+B) exp(-Ct)  and  (A/t+B) exp (-Ct) . 
A  two-term  model  consisting  of  the  sum  of  a  time  dependent 
and  time  independent  part  times  an  exponential  term  were 
required  to  fit  the  single  and  multiple  scattering  effects 
observed  in  the  data.  Two  near-surface  explosions  (E. 

Kazakh)  of  magnitude  6.1  and  one  earthquake  (N.  Argent ina- 
Chile,  559  km  depth)  of  magnitude  5.8  are  analyzed  to 
evaluate  near-source  scattering.  Five,  three-component 
North  America  RSTN  receiving  stations  are  used  to  evaluate 
the  effect  of  near-receiver  scattering. 

Coda-Q  estimates  from  the  inversions  increased  with 
frequency  in  the  approximate  range  from  800  at  1  Hz  to  1600 
at  4  Hz.  The  same  range  of  coda-Q  values  were  obtained  for 
both  of  the  models  tested;  however,  the  (A/t+B) exp (-ct) 
model  results  had  higher  variances  suggesting  poorer  fit  of 
this  model  to  the  data.  Coda-Q  estimates  are  higher  at 
stations  RSCP  and  RSNY  compared  to  the  other  stations.  The 
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coda-Q  results  from  the  GLI  analyses  are  in  agreement  with 
Q  values  determined  from  the  analysis  of  local  coda.  These 
results  suggest  that  teleseismic  coda-Q ’s  are  dominated  by 
near-receiver  scattering  effects. 

The  coda  analysis  method  outlined  in  this  paper 
represent  an  alternative  to  normal  single  scattering 
analysis  methods  performed  over  small  time  windows  early  in 
the  coda.  This  method,  which  accounts  for  multiple 
scattering  effects,  models  the  entire  coda  (500  to  600 
seconds)  where  signal  is  above  background  noise.  Good 
resilts  were  obtained  for  both  earthquake  and  explosion 
data.  While  coda  analysis  on  multiple  scattered  waves  at 
teleseismic  distances  was  the  emphasis  in  this  paper  these 
methods  should  also  work  at  local  and  regional  distances. 

INTRODUCTION 

Because  it  is  not  always  possible  to  record  data 
close  to  the  source,  coda  analysis  techniques  have  begun  to 
be  applied  to  data  at  teleseismic  distances.  At  these  large 
source-receiver  separations  multiple  scattering  effects 
play  a  dominant  role  in  coda  generation.  As  a  result  of 
this  observation  a  study  was  undertaken  to  better 
understand  how  to  deal  with  multiple  scattering  information 
in  coda  analysis.  This  is  the  first  of  three  papers  which 
address  the  issue  of  how  to  deal  with  multiple  scattering 
effects  in  coda  analysis.  Coda-Q  estimation  is  the  main 
subject  in  this  first  paper.  The  second  paper  (Soroka, 

1991) ,  referred  to  as  paper  2  in  the  text,  covers  the  issue 
of  solution  bias  due  to  model  choice,  type  of  scattering 
present  in  the  data  and  random  noise.  The  third  paper 
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(Soroka,  Stump  and  Dainty,  1991) ,  referred  to  as  paper  3  in 
the  text,  is  an  extension  of  the  typical  coda  analysis  and 
deals  with  the  conversion  of  the  A  and  B  model  parameters 
to  earth  turbidity  and  source  spectral  estimates. 

Past  work  on  seismic  coda,  the  random  energy  that 
follows  a  body  wave  arrival,  has  focused  on  data  at  local 
and  regional  distances.  This  paper  concentrates  on  coda 
following  the  direct  P-wave  arrival  with  large  source  to 
receiver  separations  in  the  range  from  66  to  98  degrees. 
Three  representative  teleseismic  waveforms  from  two 
explosions  and  one  earthquake  are  given  in  figure  1.  Shown 
are  approximately  140  seconds  of  coda  with  their 
characteristic  exponential  decay.  This  energy  travels  from 
source  to  receiver  via  many  ray  paths  and  represents 
variations  in  material  properties  along  the  way.  Because  of 
their  number  and  complexity  the  deterministic  approach  of 
locating  and  characterizing  all  acoustic  impedance 
contrasts  contributing  to  these  waveforms  is  not  feasible. 

A  statistical  approach  to  the  analysis  of  seismic  data  can 
be  applied  by  assuming  that  secondary  waves  are  generated 
at  the  numerous  random  heterogeneities  upon  incidence  of 
the  primary  wave.  The  scattered  waves  which  form  the  coda 
will,  therefore,  be  a  superposition  of  these  secondary 
waves  and  may  be  regarded  as  the  sum  of  many  independent 
small  events  (Aki,  1969). 

Past  studies  (Aki,  1969;  Aki  and  Chcuet,  1975; 

Herraiz  and  Espinosa,  1987;  and  others)  have  demonstrated 
that  coda  characteristics  for  local  earthquakes  can  be 
modeled  as  single  backscattering  processes.  Application  of 
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1.  Representative  time  series  of  the  three  events 
study.  Shown  is  the  pre-event  ambient  noise, 
direct  P-wave  arrival  and  approximately  140  sec.  of  coda. 
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the  single  backscattering  model  to  local  seismogram  codas 
has  led  to  information  concerning  the  propagation  path  and 
the  source  (Herraiz  and  Espinosa,  1987) .  Singh  and  Herrmann 
(1983)  used  the  single  backscattering  model  to  evaluate  the 
regional  characteristics  of  coda-Q  in  the  United  States. 

Aki  and  Chouet  (1975),  and  Aki  (1980a, b)  estimated  the 
turbidity  coefficient  to  obtain  information  about  the 
heterogeneities  responsible  for  coda  generation.  Aki  and 
Chouet  (1975),  Chouet  et  al.,  (1978)  and  Rautian  and 
Khalturin  (1978)  used  the  coda  to  obtain  estimates  of  the 
source . 

Frankel  and  Clayton  (1984)  used  finite  difference 
modeling  techniques  and  found  that  the  apparent  attenuation 
due  to  scattering  was  greatest  when  the  scatter er  size  was 
comparable  to  the  seismic  wavelength.  Smaller  scale 
features  were  found  to  have  an  important  but  lesser  effect 
on  coda  generation.  An  estimate  of  the  seismic  wavelength 
can  be  computed  from  the  dominant  frequency  of  the  source 
wavelet  and  the  velocity  of  the  material  in  which  the  wave 
is  propagating.  For  the  data  of  this  study,  the  estimated 
wavelengths  range  from  1000  to  5000  meters  for  frequencies 
from  1  to  3  Hz  and  velocities  from  3  to  5  km/sec.  Low 
velocities  are  used  because  past  coda  work  suggests  shear 
wave  scattering  is  responsible  for  coda  generation  (Herraiz 
and  Espinosa,  1987) . 

Strong  evidence  for  low  apparent  phase  velocities 
(3. 5-4. 5  km/sec)  in  teleseismic  P-wave  codas  for  deep 
earthquakes  and  explosions  were  observed  by  Dainty  (1985) 
at  NORSAR  and  NORESS.  These  phase  velocities  were 
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determined  from  wavenumber  analysis  of  vertical  component 
seismograms  and  interpreted  as  teleseismic  P  to  Lg  (trapped 
shear  wave)  scattering  near  the  receiver.  For  the 
explosions  Dainty  (1985)  also  found  an  equal  amount  of  high 
phase  velocities  that  he  interpreted  as  Lg  to  teleseismic  P 
scattering  near  the  source.  The  data  analyzed  in  this  study 
show  a  significant  amount  of  energy  on  the  radial  and 
transverse  components  suggesting  P  to  Lg  scattering  near 
the  receiver. 

The  teleseismic  coda  analysis  used  in  this  paper  is 
unique  in  a  number  of  respects.  The  modeling  includes  both 
single  and  multiple  scattering  effects.  Generalized  Linear 
Inversion  methods  are  applied  so  that  the  entire  coda  where 
signal  is  sufficiently  above  noise  can  be  used.  Inversion 
techniques  are  employed  to  handle  large  quantities  of  data 
to  improve  the  resolution  of  model  parameters  as  a  function 
of  frequency.  The  larger  data  set  reduces  the  influence  of 
localized  data  problems  that  might  bias  the  final 
interpretations . 

A  review  of  pertinent  past  coda  analysis  work  will  be 
presented  first.  This  is  followed  by  a  description  of  the 
methods  used  in  this  study.  The  results  are  described  last 
and  include  comparisons  to  other  independent  measures  of 
similar  properties. 


BACKGROUND 

Aki  and  Chouet  (1975)  successfully  used  single 
backscattering  theory  to  model  the  codas  of  local  seismic 
events.  The  wide  applicability  of  the  single  backscattering 

coda  model  has  been  demonstrated  by  its  use  in  local  and 
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regional  coda  analysis  (Herraiz  and  Espinosa,  1987) .  These 
studies  include  the  use  of  codas  to  deterinine  attenuation 
or  coda-Q  (Aki  and  Chouet  (1975),  Rovelli  (1982),  Singh  and 
Herrmann  (1983)),  the  source  spectrum  (Aki  and  Chouet 
(1975),  Chouet  et  al.,  (1978),  Rautian  and  Khalturin 
(1978))  and  turbidity  (Aki  and  Chouet  (1975),  Aki 
( 1980a, b).  Dainty  et  al.,  (1987)).  The  potential  for  using 
teleseismic  coda  in  similar  fashion  was  demonstrated  by  Aki 
(1982).  The  codas  were  found  to  be  insensitive  to  the  path 
and  to  have  similar  amplitude  and  spectral  characteristics 
at  different  stations  for  the  same  source.  The  teleseismic 
codas  were  modeled  as  a  superposition  of  many  random 
independent  single  backscattered  events  similar  to  the 
methods  used  in  local  and  regional  coda  analysis. 

The  single  backscattering  model  describes  the  time 
(t)  and  frequency  (f)  variation  of  the  power  spectrum  of 
the  coda  (P(f,t)),  (Herraiz  and  Espinosa,  1987).  The  model 
parameters  are  the  source  term  (S(f)),  coda-Q  (Q(f))  which 
includes  intrinsic  attenuation  and  scattering  attenuation 
effects  and  turbidity  (T(f)).  Turbidity  is  a  measure  of  the 
medium's  ability  to  initiate  scattering.  It  represents  the 
wave  path  heterogeneity  and  is  expressed  as  a  cross-section 
per  unit  volume  with  units  of  1/km.  The  mathemat- ^  cal 
expression  for  the  power  spectrum  of  single  backscattered 
body  waves  (Gao  et  al.,  (1983b))  is: 

2  T(f)  s(f)  -2^ft 

Pl(f,t)  =  - - -  exp( - )  (1) 

V  t2  Q(f) 

where:  P(f,t)  =  power  spectrum 

S(f)  =  source  term  =  r^  PQ(f)  exp(47rfr/VQ) 
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T(f)  =  turbidity 

Q(f)  =  coda-Q 

V  =  average  crustal  velocity 

r  =  reference  distance 

The  single  backscattering  coda  model  is  not 
applicable  to  all  seismic  data.  In  some  cases  evidence  for 
multiple  scattering  effects  have  been  presented  (Richards 
and  Menke,  1983) .  In  an  attempt  to  better  account  for  these 
multiple  scattering  effects  the  energy  flux  model  has  been 
applied  (Frankel  and  Wennerberg  (1987)  and  Langston 
(1989)).  Gao  et  al.,  (1983b)  extended  the  single 
backscattering  model  to  multiple  scattering  and  provided 
the  mathematical  expressions  of  higher  order  scattering 
models.  They  showed  that  single  scattering  effects  dominate 
at  early  times  in  the  coda  and  multiple  scattering  with  its 
slower  decay  become  important  at  longer  times.  The  slower 
decay  rate  of  multiple  scattering  is  a  result  of  the  time 
delay  of  primary  wave  energy  by  scattering  interactions. 

The  slower  decay  rate  associated  with  the  multiple 
scattering  process  can  explain  the  long  coda  durations 
which  are  characteristic  of  teleseismic  observations 
(figure  1).  The  mathematical  expressions  for  second  and 
third  order  multiple  backscattering  (Gao  et  al.,  (1983b)) 
are: 


-27rft 

P2(f,t)  = 

2.46  srf)  V 
t 

Tff)  2 

exp( 

- ) 

Q(f) 

-27rft 

(2) 

P3(f.t)  = 

0.716  S(f) 

T(f)  2 

V  exp( 

- ) 

Q{f) 

(3) 

A  simple  model  which  consists  of  a  time  dependent  and 
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time  independent  term  of  the  form  (A/t^+B) exp(-Ct)  is  seen 
to  correspond  to  first  and  third  order  scattering  in 
equations  1  and  3  above.  Attempts  were  made  to  include 
second  order  scattering  effects  but,  due  to  the  extra 
degree  of  freedom  the  inversions  were  unstable  on  all  but 
noise-free  synthetics.  First  plus  third  order  scattering 
was  found  to  be  a  good  average  model  to  fit  to  the  data. 

The  results  of  these  model  tests  are  described  in  more 
detail  in  paper  2.  Note  that  another  simple  model 
consisting  of  a  time  dependent  and  time  independent  term  is 
of  the  form  (A/t+B) exp(-Ct) .  This  alternative  model  also 
describes  the  characteristics  of  the  data. 

Representative  teleseismic  coda  spectral  decay  curves 
at  a  single  frequency  (2  Hz)  are  shown  in  figure  2.  The 
data  are  seen  to  decay  rapidly  at  early  times  and  more 
slowly  at  later  times.  Superimposed  on  the  data  is  the 
least-square  best  fit  of  the  (A/t^+B) exp(-Ct)  model  in  the 
upper  display  and  the  (A/t+B) exp (-Ct)  model  in  the  lower 
display.  The  time  independent  multiple  scattering  teirm  fits 
the  data  well  at  later  times  but  has  poor  predictions  at 
early  times.  The  time  dependent  single  scattering  model 
curve  fits  the  data  well  at  early  times  and  diverges  at 
later  times.  Attempts  to  use  the  time  dependent  single 
scattering  term  to  describe  the  entire  da  za  set  failed 
because  time-varying  coda-Q's  were  required  to  effectively 
model  the  data  (Soroka  et  al.,  (1985)).  These  results 
indicate  that  a  two  term  model  which  includes  a  time 
dependent  single  and  time  independent  multiple  scattering 
term  is  necessary  to  describe  coda  when  multiple  scattering 
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Figure  2.  Representative  (A/t^)  and  (A/t)  model  fits  to 
the  data.  The  single  and  multiple  models  are  plotted 
separately  to  illustrate  that  neither  of  the  models  alone 
can  describe  the  data  at  all  times.  Analysis  is  for  the 
vertical  component  at  2  Hz  of  event  2  at  station  RSON. 
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effects  are  present.  This  two-term  model  appears  more 
appropriate  than  the  single  scattering  model  alone  for 
describing  teleseismic  coda  and  probably  seismic  coda  in 
general  when  long  time  windows  which  include  multiple 
scattered  energy  are  used  (paper  2) . 

In  summary  the  two  models  which  adequately  fit  the 
time  varying  characteristics  observed  in  the  data  are: 

A 

(  -  +  B  )  exp(-Ct)  (4) 

t2 

and: 

A 

(  -  +  B  )  exp(-Ct)  (5) 

t 

These  models  will  be  referred  to  in  the  text  as  (A/t^)  and 
(A/t)  respectively.  The  (A/t^)  model  is  shown  above  to  be 
the  form  of  single  and  triple  scattered  body  waves  as 
derived  by  Gao  et  al.,  (1983b). 

The  (A/t)  model  is  the  form  of  two  different  physical 
processes,  one  more  and  one  less  applicable  to  the 
teleseismic  coda  analysis  problem.  The  time  dependent  and 
time  independent  terms  in  the  (A/t)  model  correspoj-.d  to 
first  and  second  order  surface  wave  scattering  (Gao  et  al., 
1983a) .  Evidence  presented  later  in  this  paper  favors  body 
wave  over  surface  wave  scattering  in  teleseismic  coda 
generation.  The  time  dependent  term  in  the  (A/t)  model  is 
also  the  correct  form  of  single  scattered  body  waves  that 
result  from  a  plane  wave  encountering  a  layer  of  randomly 
distributed  scatterers.  The  time  independent  term  could 
represent  trapped  P-to-Lg  waves  that  are  generated  when  the 
plane  wave  hits  rhe  near  surface  layer  under  the  receiver. 
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The  plane  wavs  model  is  discussed  in  more  detail  in  paper 
J.  For  the  purposes  of  coda-Q  estimation  it  is  only 
important  to  note  that  the  exponential  term  in  both  the 
(A/t^)  and  A/t)  models  (equations  4  and  5)  are  the  same. 
Coda-Q  estimates  will  be  determined  with  both  models  and 
the  results  compared  in  the  detailed  result  section. 

PROCEDURES 

The  coda  analysis  method  used  in  this  study  involves 
fitting  a  model  to  time  varying  spectral  estimates  of 
teleseismic  coda  and  then  computing  earth  and  source 
properties  from  the  model  parameters.  Both  of  the  proposed 
models  (equations  4  and  5)  are  fit  to  the  data  using 
generalized  linear  inversion  (GLI)  techniques.  The  analysis 
procedure  is  divided  into  four  steps: 

1)  TIME  VARYING  SPECTRAL  CHARACTERIZATION 

2)  CODA  MODEL  VERIFICATION 

3)  INVERSION 

4)  INTERPRETATION 

Each  of  these  four  topics  will  be  discussed  in  turn. 
First  however,  the  selection  of  an  appropriate  data  set  for 
analysis  will  be  considered.  The  data  chosen  possessed  uhe 
following  characteristics; 

1)  Well  developed  codas  following  the  direct  P-wave 

2)  Events  recorded  at  teleseismic  distances 

3)  Coda  uncorrupted  by  deterministic  arrivals 

4)  Digitally  recorded  data  with  a  sample  rate  that  allows 

for  high  frequency  (1-10  Hz)  analysis 
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Data  from  the  Regional  Seisr..ic  Test  Network  (RSTN)  , 
(Breding,  1983)  seismic  database  (table  1)  met  these 

Table  1.  — Event  and  station  information 


EVENT  LOCATIONS 


EVENT 

TYPE 

LOCATION 

1 

EXPLOSION 

EASTERN  KAZAKH,  U.S.S.R. 

2 

EXPLOSION 

EASTERN  KAZAKH,  U.S.S.R. 

3 

EARTHQUAKE 

N.  ARGENTINA-CHILE 

STATION  LOCATIONS 

STATION 

LOCATION 

RSCP 

CUMBERLAND  PLATEAU,  U.S.A. 

RSNY 

NEW  YORK,  U.S.A. 

RSSD 

SOUTH  DAKOTA,  U.S.A. 

RSON 

ONTARIO,  CANADA 

RSNT 

NORTHWEST  TERRITORY,  CANADA 

SOURCE-TO  RECEIVER  DISTANCES  IN  DEGREES 


EVENT  MAGNITUDE  DEPTH  STATION 


RSNT 

RSON 

RSNY 

RSSD 

RSCP 

1 

6.1 

NEAR-SURFACE 

68 

79 

83 

86 

94 

2 

6.1 

NEAR-SURFACE 

68 

79 

83 

86 

94 

3 

5.8 

559  Km 

98 

82 

72 

80 

66 

requirements.  Two  near-surface  explosions  (mj3=6.1)  and  one 
deep  earthquake  (559  km  depth,  mj3=5.8)  were  selected  for 
analysis.  Shallow  earthquakes  were  not  chosen  because  the 
pP  (surface  reflection)  occurs  too  soon  after  the  direct 
P-wave  arrival  and  corrupts  the  coda  decay.  Two  explosions 
and  one  deep  earthquake  were  selected  for  the  purpose  of 
studying  source  differences  and  near-source  versus 
near-receiver  scattering  effects.  Records  from  five  North 
American  stations  with  source-to-receiver  distances  that 
range  from  60  to  98  degrees  were  used  to  study  station 
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differences.  The  RSTN  data  were  digitally  recorded  on  three 
components  (vertical,  north  and  east)  at  a  sample  rate  of 
40  samples  per  second  (Breding,  1983) .  Seismograms  for  the 
three  events  at  all  five  North  America  RSTN  stations  (table 
1)  were  analysed.  The  north  and  east  components  are 
rotated  to  form  radial  and  transverse  components  (to  the 
source)  prior  to  analysis. 

TIME  VARYING  SPECTRAL  CHARACTERIZATION 
The  power  spectra  were  computed  using  the  Fast 
Fourier  Transform  (FFT)  method  on  small,  3.75  second  time 
windows  (150  samples).  To  reduce  edge  effects,  a  Hanning 
window  was  applied  to  the  selected  data  with  zero  padding 
prior  to  performing  the  FFT.  Figure  3  is  a  representative 
3-D  display  of  the  power  spectral  decay  characteristics 
with  frequency  and  time.  Plotted  is  the  signal-to-noise 
ratio  computed  from  the  data  and  average  pre-event  ambient 
noise  spectral  estimates.  All  3-D  disp3ays  of  the  time 
varying  spectra  have  certain  features  in  common.  Maximum 
spectral  power  occurs  at  the  direct  arrival  between  1 
and  3  Hertz.  Spectral  amplitude  decreases  with  time  at  all 
frequencies.  The  signal  falls  below  background  noise  sooner 
at  high  and  low  frequencies  relative  to  the  spectral 
maximum . 


CODA  MODEL  SELECTION 

Selection  of  the  appropriate  model  for  the 
teleseismic  coda  is  not  generally  clear.  The  approach  by 
some  authors  is  to  use  the  single  backscattering  model  for 
the  sake  of  consistency  with  previous  work.  While  this 
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Figure  3.  Time  varying  power  spectra  of  Event  2.  Shown  is 
the  signal-to-noise  ratio  in  db  for  the  vertical  component 
at  station  RSNT. 
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approach  is  the  simplest  it  can  lead  to  biased  solutions 
which  may  not  reflect  true  earth  properties.  There  are 
three  main  issues  that  need  to  be  addressed  when 
determining  which  model  is  most  appropriate  for  the  coda 
under  analysis. 

1-  FORWARD  -VERSUS-  BACK-SCATTERING 

2-  SINGLE  -VERSUS-  MULTIPLE  SCATTERING 

3-  SURFACE  -VERSUS-  BODY  VJAVE  SCATTERING 

FORWARD  -VERSUS-  BACK-SCATTERING 

Menke  and  Chen  (1984)  and  Richards  and  Menke  (1983) 
have  proposed  criteria  for  distinguishing  between  forward 
and  backscattering  effects  in  coda.  They  include: 

For  forward  scattered  wo'-.s: 

a)  the  coda  has  relatively  higher  frequencies  than  the 

initial  wave.  (Richards  and  Menke,  1983) 

b)  the  coda  power  spectra  deca^  with  time  is  faster 

for  lower  frequencies  than  highf'r  frequencies. 

(Menke  and  Chen,  1984;  Frankel  and  Clayton,  1984) 

For  backscattered  waves; 

a)  at  all  frequencies  the  coda  spectral  power  is  lower 

than  the  initial  wave.  (Richards  and  Menke,  1983) 

b)  the  coda  power  spectra  decay  rate  is  faster  for  high 

frequencies  then  low  frequencies. 

(Menke  and  Chen,  1984;  Frankel  and  Clayton,  1984) 

c)  estimates  of  apparent  Q  made  from  coda  of  backscattered 

waves  increase  with  frequency,  those  from  forward 
scattered  waves  do  not  (Richards  and  Menke,  1983). 
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spectral  ratios  (Pizarenko,  1S‘'0)  were  used  to 
characterize  the  changes  in  relative  amplitude  of  the  coda 
spectrum  with  time  as  a  function  of  frequency.  Figure  4  is 
a  representative  spectral  ratio  plot  (Pj^  /  Pi-i  ,  where  P^ 
is  the  spectral  power  at  a  specific  time)  and  shows  that 
power  decreases  with  time  at  all  frequencies.  There  is  no 
appreciable  increase  in  power  at  high  frequencies  after  the 
initial  wave  that  cannot  be  attributed  to  high  variance  in 
the  spectral  estimates.  The  coda-Q  solutions,  reported 
later,  are  also  found  to  increase  with  frequency  in 
agreement  with  the  backscattering  model  criteria.  These 
observations  lead  to  the  conclusion  that  teleseismic  coda, 
for  the  three  events  under  analysis,  are  dominated  by  the 
backscattering  processes . 

SINGLE  -VERSUS-  MULTIPLE  SCATTERING 

It  was  shown  earlier  for  both  the  (A/t^)  and  (A/t) 
models  in  figure  2  that  neither  the  single  or  multiple 
model  alone  described  the  characteristics  of  the  data.  The 
time  dependent  single  scattering  curve  fit  the  coda  at 
early  times  and  the  time  independent  multiple  scattering 
curve  fit  the  data  at  later  times.  Figure  5  shows  that  goo:, 
agreement  is  obtained  between  the  data  and  both  models  when 
the  single  and  multiple  terms  are  added  prior  to  being 
overlain  on  the  data.  The  (A/t^)  and  (A/t)  models  are 
plotted  on  top  of  the  time  varying  spectral  data  and  can  be 
seen  to  follow  the  trend  of  the  data  at  all  times.  The  only 
exception  might  be  at  very  early  times  where  the  (A/t) 
model  with  its  gentler  decay  does  not  follow  the  data  as 
well  as  the  (A/t^)  model.  This  small  difference  between  the 
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Figure  5.  Representative  (A/t^)  and  (A/t)  final  model  fits 
(lines)  to  the  data  (points) .  The  models  plotted  are  the 
single  plus  multiple  models  discussed  in  the  text.  A  five 
point  smoothing  function  is  applied  to  reduce  the  high 
frequency  noise  in  the  data.  Analysis  is  for  the  vertical 
component  at  2  Hz  of  Event  2  at  station  RSON. 
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two  models  at  early  times  in  the  coda  may  not  be 
resolvable.  On  the  other  hand,  results  from  paper  2 
Indicate  that  this  feature  may  make  one  model  less  stable 
than  the  other  during  the  inversions.  Based  on  these  model 
fits  to  the  data  a  model  consisting  of  a  time  dependent 
single  scattering  term  and  a  time  independent  multiple 
scattering  term  appears  necessary  to  model  coda  with 
multiple  energy  present. 

Attempts  were  made  to  perform  the  inversions  with 
additional  orders  of  scattering  (or  additional  time 
dependent  terms)  in  the  hope  that  this  would  improve  fits 
and  reduce  variances.  The  details  of  these  tests  are 
included  in  paper  2.  These  additional  free  parameters  were 
found  unstable  on  all  but  noise  free  synthetic  data.  The 
inversion  technique  is  only  able  to  resolve  the  steeply 
decaying  early  time  coda  from  the  gently  sloping  later  time 
coda.  The  inversion  was  also  found  to  be  very  unstable  when 
the  input  data  contained  orders  of  scattering  significantly 
different  from  the  model  being  fit. 

In  summary,  the  decision  to  use  the  (A/t^)  and  (A/t) 
models  is  based  on: 

1)  the  good  agreement  between  data  and  models,  the  two-term 

models  were  found  to  be  good  average  models  to  use 

2)  the  necessity  of  keeping  the  model  simple  for  use  with 

the  inversion  techniques 

3)  the  inability  of  the  inversion  tests  using  additional 

levels  of  scattering  to  improve  the  fits 

SURFACE  -VERSUS-  BODY  WAVE  SCATTERING 

Distinguishing  surface  from  body  wave  scattering  is 
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not  something  that  is  easily  tested.  The  evidence  in  the 
literature  appears  to  strongly  favor  body  wave  scattering 
as  the  mechanism  for  coda  generation.  Herraiz  and  Espinosa 
(1987)  use  the  following  three  arguments  in  support  of  body 
wave  scattering; 

1)  Coda  waves  corresponding  to  surface  and  a  deep 

borehole  ('3.5  km)  site  in  Japan  share  the  same 

features  (Sato,  1978) . 

2)  Coda  and  S-waves  have  similar  site  effects  for  a 

wide  frequency  range  1  to  25  Hz  (Tsujiura,  1978) . 

3)  Q  for  shear  waves  and  coda-Q  are  similar  over  che 

frequency  range  1  to  25  Hz  (Aki,  1980a, b) . 

Based  on  the  Gao  et  al.,  (1983a  and  b)  theory  the 
difference  between  the  body  and  surface  wave  models  is  in 
the  time  dependence  of  the  single  scattering  term.  For  body 
waves  it's  A/t^  and  for  surface  waves  A/t.  Both  of  these 
models  fit  the  time  varying  spectra  of  teleseismic  coda 
(figure  5)  and  produce  stable  solutions.  Based  on  the  above 
three  pieces  of  evidence  surface  waves  can  be  ruled  out  as 
a  major  contributor  in  teleseismic  coda  generation. 

Surface  wave  coda-Q  estimates  would  be  identical  to 
those  deteirmined  by  fitting  the  (A/t)  model  which  also 
describes  a  plane  wave  hitting  a  layer  of  random  scatterers 
and  is  included  in  this  study.  The  important  difference 
between  the  surface  wave  and  plane  wave  models  is  with  the 
definition  of  the  A  and  B  model  parameters  which  are 
different  for  each  model.  The  exponential  (C)  term  from 
which  coda-Q  is  derived  would  be  the  same  in  each  case 
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(paper  3).  The  coda-Q  results  from  fitting  the  (A/t^)  and 
(A/t)  models  will  be  presented  in  the  uetailed  results 
section  for  comparison.  The  conversion  of  the  A  and  B  model 
parameters  to  earth  and  source  property  et.timates  is  the 
subject  of  paper  3. 

GENERALIZED  LINEAR  INVERSION 

The  scattering  model  parameters  are  determined  from 
the  time  varying  spectral  data  using  the  GLI  technique  (Aki 
and  Richards,  1980  and  Henke,  1984) .  Frequencies  as  high  as 
10  Hz  are  analyzed  for  each  event,  station  and  component. 

A  signal-to-noise  ratio  of  10  db  and  a  minimum  of  10  data 
values  were  required  for  the  inversion.  The  first  of  the 
two  parametrized  models  used  in  the  inversion  is: 

A 

P  =  (  +  B  )  EXP(  -C  t  )  (6) 

t2 

Comparison  of  this  model  to  the  Gao  et  al.,  (1983b) 
theory,  equations  (1)  and  (3)  give  the  following 
relationship  between  the  model  parameters  A  and  B  and  Earth 
and  source  properties  described  in  the  introduction; 

A=2ST/V,  B  =  0.716  S  T^  V  ,  C  =  27rf/Q  (7) 

The  emphasis  of  this  paper  is  on  the  estimation  of 
coda-Q  from  teleseismic  data  where  multiple  energy  is  a  ' 

dominant  contributor  to  the  coda.  Equation  8  below  gives 
the  alternative  model  which  also  describes  the 
characteristics  of  the  data.  The  model  parameter  C  in 
equation  (8)  from  which  coda-Q  is  computed  is  the  same  as 
for  equation  (6)  (paper  3).  Coda-Q  is  therefore  obtained 
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from  C  in  the  same  manner  described  above  (C=27rf/Q)  . 

P  *  (  A/t  +  B  )  exp(-C  t)  (8) 

The  A  and  B  model  parameters  are  not  the  same  for  the 
(A/t^)  and  (A/t)  models.  A  detailed  study  of  turbidity  and 
source  spectral  estimation  from  the  A  and  B  model 
parameters  is  the  subject,  of  paper  3.  Using  the  (A/t^) 
modfcl,  which  appears  to  fit  the  data  better,  a  limited  test 
was  performed  in  which  turbidity  is  computed  from  the  ratio 
of  B/A  and  assuming  a  material  velocity  of  3.5  km/sec.  The 
source  term  is  obtained  from  A.  The  results  of  this  limited 
test,  shown  in  the  result  section,  suggest  that  the  simple 
two-term  representation  of  the  complicated  coda  scattering 
problem  could  provide  a  means  of  estimating  turbidity  and 
source  spectral  properties.  The  reader  is  referred  to  paper 
3  for  more  details  on  this  work. 

ERROR  ANALYSIS 

An  extensive  study  was  undertaken  to  better 
understand  the  significance  of  the  coda  analysis  results 
when  multiple  scattered  energy  is  present.  The  details  of 
this  study  are  the  subject  of  paper  2.  A  brief  summary  of 
the  findings  are  presented  below  and  the  reader  is  referred 
to  paper  2  for  more  detail. 

Using  synthetics,  paper  2  examines  the  effects  of 
multiple  scattered  energy,  random  noise  and  model  parameter 
magnitude  on  inversion  stability  and  solution  accuracy.  The 
popular  models  used  in  coda  analysis  are  compared  to  better 
understand  their  similarities,  differences  and  limitations. 
Also  covered  is  the  issue  of  solution  bias  which  results 
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from  using  simple  two-term  models  to  represent  different 
types  of  multiple  scattered  energy  in  the  coda. 

A  Monte  Carlo  error  analysis  was  run  first  to  test 
the  effect  of  random  noise  on  the  inversion  results.  If 
random  noise  effects  severely  corrupted  the  inversion 
solutions  there  was  no  point  in  proceeding.  The  test 
consisted  of  fitting  the  (A/t^)  model  100  times  to  a  known 
input  model  while  adding  different  sets  of  random  noise  of 
a  constant  level  to  the  input  data.  The  noise  added  to  the 
data  was  produced  by  generating  a  uniform  distribution  of 
random  numbers.  The  four  levels  of  noise  tested  are 
approximately  1,  10,  20  and  40  percent  of  signal  maximum. 
The  coda-Q,  turbidity  and  source  term  results  from  the  100 
runs  were  used  to  calculate  a  mean  and  standard  deviation 
for  each  of  these  parameters.  Figure  6  illustrates  the 
solution  uncertainty  (%  error=standard  deviation  in  percent 
of  mean)  for  larger  values  of  coda-Q,  turbidity,  and  source 
term  as  a  function  of  random  noise. 

The  difference  betv/een  the  Monte  Carlo  mean  values 
and  the  known  model  is  a  measure  of  the  bias  expected  in 
the  inversions  due  to  random  noise.  Table  2  gives  the  bias 
for  each  of  the  parameters  and  shows  (similar  to  figure  6) 
that  coda-Q  and  turbidity  (T)  can  be  computed  with  greater 
certainty  than  the  source  term  (S) . 

A  similar  Monte  Carlo  test  with  low  values  of  Q, 
turbidity  and  source  term  suggested  that  the  results  in 
figure  6  and  table  2  are  probably  a  worst  oase  situation. 
The  level  of  noise  in  the  observational  data  varies  from 
station  to  station  and  with  frequency  but  was  determined  to 
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be  within  the  range  indicated  on  figure  6  (7-14%) .  These 
random  noise  tests  indicate  that  coda-Q  can  be  determined 
more  accurately  then  the  other  two  parameters  (T  and  S)  and 
Q  is  relatively  insensitive  to  even  high  levels  of  random 
noise.  Turbidity  can  be  determined  more  accurately  then  the 
source  term  and  should  produce  relatively  stable  estimates. 
The  source  term  should  have  the  largest  degree  of 
uncertainty. 

TABLE  2. — Solution  bias  due  to  random  noise  (percent 
difference  between  the  Monte  Carlo  analysis  mean  and  the 

known  solution) 


MODEL  PARAMETER 

1% 

RANDOM  NOISE 
10% 

LEVELS 

20% 

40% 

Q 

0.0 

0.L> 

3. 

10. 

T 

0.0 

3. 

10. 

115. 

S 

0.1 

30. 

200. 

2000. 

The  Monte  Carlo  analysis  indicated  that  the  effects 
of  random  noise  were  tolerable  and  that  stable  coda 
analysis  estimates  could  be  detennined.  To  better 
understand  the  applicability  of  the  simple  two-term  models 
a  comparison  was  made  between  the  (A/t^) ,  (A/t)  and  energy 
flux  models.  The  results  of  this  comparison  showed  them  to 
produce  equivalent  estimates  of  coda-Q.  However,  when  the 
single  scattering  model  was  fit  to  data  with  single  plus 
multiple  scattered  energy,  the  coda-Q  solutions  were 
strongly  biased  by  as  much  as  280%.  A  similar  situation  was 
observed  when  the  (A/t^)  model  was  fit  to  data  with  only 
single  scattered  energy.  The  results  of  these  tests  tell  us 
that  poor  coda-Q  estimates  are  obtained  when  the  model  to 
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be  fit  to  the  data  differs  significantly  from  the 
scattering  present  in  the  data.  The  two  term  models  (A/t^) 
and  (A/t)  were  found  to  be  good  average  models  of  up  to 
fourth  order  multiple  energy.  The  (A/t^)  model  was  observed 
to  described  steep  decays  at  early  times  better  than  the 
(A/t)  model. 

A  comparison  was  also  made  between  the  A  and  B  model 
parameters  from  the  (A/t^)  and  (A/t)  models.  These  results 
showed  large  differences  which  suggests  that  model 
selection  may  be  critical  if  useful  information  is  to  be 
extracted  from  the  A  and  B  model  parameters.  This  is  the 
subject  of  paper  3  to  which  the  reader  is  referred  for  more 
detail . 


INTERPRETATION  OF  MODELING  RESULTS 
The  (A/t2)  and  (A/t)  scattering  models  were  applied 
to  all  RSTN  data  from  the  two  explosions  and  one  earthquake 
(table  1).  Figures  7  and  8  are  representative  (A/t^)  and 
(A/t)  model  solutions  of  coda-Q  versus  frequency  for  the 
three  events  recorded  at  station  RSON.  Vertical,  radial  and 
transverse  solutions  are  superimposed  on  each  plot.  In  each 
case  the  results  are  observed  to  follow  a  linear  increase 
of  coda-Q  with  frequency.  The  (A/t^)  and  (A/t)  model 
results  are  similar  but  not  equal.  The  increased  scatter  in 
coda-Q  at  low  and  high  frequencies  is  due  to  low 
signal-to-noise  ratios.  Consistent  coda-Q  estimates  from 
the  three  components  support  the  conclusion  that  each 
component  is  measuring  similar  earth  properties.  The 
agreement  in  coda-Q  estimates  between  the  deep  earthquake 
and  the  two  shallow  explosions  suggests  that  near-receiver 
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Figure  7.  (A/t^)  model  coda-Q  versus  frecjuency  results  for 

Events  1,  2  and  3  at  station  RSON.  The  slope  and 
Y-intercept  values  for  the  least-square  best  fit  lines  are 
shown  on  the  figures.  Symbols  represent  the  different 
components:  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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Figure  8.  (A/t)  model  coda-Q  versus  frequency  results  for 

Events  1,  2  and  3  at  station  RSON.  The  slope  and 
Y-intercept  values  for  the  least-square  best  fit  lines  are 
shown  on  the  figures.  Symbols  represent  the  different 
components;  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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scattering  dominates  over  near-source  scattering  in  these 
codas.  The  implications  of  this  observation  are  that  coda 
analysis  reflect  earth  properties  near  the  receiver. 

Figure  9  is  a  representative  display  of  both  the  A 
and  B  model  parameters,  from  equation  4,  as  a  function  of 
frequency.  These  displays  do  not  tell  us  about  earth  and 
source  properties  but  suggest  that  some  systematic  changes 
are  occurring  in  the  data  and  are  being  measured  by  these 
parameters.  To  extract  earth  and  source  properties  from  A 
and  B  it  is  necessary  to  interpret  A  and  B  in  terms  of  a 
scattering  theory.  Equations  1  and  3  provide  a  means  of 
converting  A  and  B  into  turbidity  (T)  and  source  spectra 
(S)  by  equating  the  A  parameter  to  the  single  scattering 
amplitude  term  (in  front  of  the  exponential)  and  B  to  the 
triple  scattering  amplitude  term.  The  extraction  of  earth 
and  source  information  from  the  A  and  B  model  parameters  is 
a  potential  new  source  of  coda  information.  The  issues  of 
proper  model,  model  assumptions  and  estimate  reliability  is 
discussed  in  detail  in  paper  3. 

Because  the  (A/t^)  model  allows  for  the  separation  of 
the  turbidity  and  source  term,  a  small  test  was  performed 
with  a  (A/t^)  model  to  determine  if  additional  work  in  this 
area  was  warranted.  In  this  test  turbidity  and  source 
spectra  are  computed  from  A  and  B  for  a  representative  data 
set.  The  resulting  estimates  are  compared  with  other 
independent  measures  of  these  properties  and  evaluated  in 
terms  of  general  appearance. 

Figure  10  is  a  representative  plot  of  turbidity 
(T(f))  for  all  three  events  at  station  RSNT  with  the  three 


no 


^ '  _  Model  parameter  A 

8.  - 


7.  H  . 


0.10  0.25  0.63  1.6  4.0  10.0 


FREQUENCY  (Hz) 


FREQUENCY  (Hz) 


^  ®  model  parameter  solutions 

from  the  (A/t^+B) exp(-Ct)  model. 


to  0 


EVENT  1 


y  0 

>» 

0  10 

5 

1 

L4  _ 1 

S 

X 

0  0,010  1 

H  I 

0  0010 

Slop*  s  *  <001  ±  2<*% 

/s  ft  A 

tnt*re*pt  *  ♦  .023  ±  6% 

20  40  60 

FREQUENCY  (HZ) 


EVENT  2 


EVENTS 


Figure  10.  (A/t^)  model  turbidity  versus  frequency  results 

for  Events  1,  2  and  3  at  station  RSSD.  The  least-square 
best  fit  line  is  plotted.  Symbols  represent  the  three 
components:  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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components  superimposed.  Turbidity,  which  is  a  measure  of  a 
medium's  ability  to  initiate  scattering,  appears  to  be 
relatively  constant  with  frequency.  The  agreement  among  the 
different  components  suggests  they  are  again  representative 
of  similar  earth  properties.  The  scatter  in  the  turbidity 
results  tends  to  be  higher  then  with  the  Q  results  as 
predicted  by  the  Monte  Carlo  error  analysis.  The  values 
fall  consistently  between  .01  and  .1  which  is  in  agreement 
with  other  published  work  (Aki,  1973,  Dainty,  1988). 

Because  of  the  high  variance  the  gentle  slopes  of  the  best 
fit  least-square  lines  are  probably  not  significant. 

Figure  11  is  a  representative  display  of  the  source 
term  result  for  the  three  events  at  station  RSON  (three 
components  superimposed) .  The  vertical  scale  is  spectral 
power  in  db  down  from  maximum.  The  maximum  spectral  value 
used  to  compute  the  db  spectrum  is  constant  for  a 
pa:  -icular  source,  receiver  and  component  but  varies 
otherwise  to  make  the  results  more  comparable.  The  expected 
high  variance  in  the  results  has  been  suppressed  because 
log  scales  are  used  for  both  axes.  Despite  the  noise 
problems  the  source  spectrums  computed  with  the  simplified 
theory  have  the  correct  general  appearance  above  1  Hz  where 
the  signal-to-noise  ratio  is  better. 

Since  pre-event  noise  was  subtracted  from  the  data 
prior  to  analysis,  low  signal-to-noise  ratio  values  will 
strongly  effect  the  shape  of  the  final  source  spectrums. 
Pre-event  noise  was  relatively  flat  above  1  Hz  but  rose  to 
high  values  below  1  Hz.  As  a  result  the  source  term 
determinations  below  1  Hz  are  strongly  affected  while  those 
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Figure  ll.  (A/t^)  model  source  term  solutions  for  Events 
1,  2  and  3  at  station  RSON.  Plot  is  spectral  power  in  db 
down  from  maximum  with  the  instrument  response  removed. 
Symbols  represent  the  three  components:  vertical  (+) , 
radial  (x)  and  transverse  (o) . 
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above  1  Hz  are  not. 

Based  on  the  encouraging  results  from  this  small  test 
a  more  extensive  study  was  undertaken  in  which  both  the 
inversion  results  from  the  (A/t^)  and  (A/t)  models  are 
used.  Details  of  this  work  can  be  found  in  paper  3  and  will 
not  be  discussed  any  further  here. 

DETAILED  RESULT  COMPARISONS 

Three  types  of  composite  displays  are  provided  for 
the  coda-Q  estimates.  The  first  is  a  plot  by  common 
component  (vertical,  radial  and  transverse) .  The  second  is 
by  common  event  (1,  2  and  3)  and  finally  by  common 
receiving  station  (RSON,  RSNT,  RSSD,  RSCP  and  RSNY) .  The 
(A/t^)  and  (A/t)  model  results  for  each  of  the  3  composite 
displays  will  be  shown  together  for  easier  comparison. 

The  coda-Q  versus  frequency  results  are  plotted  by 
common  component  for  the  (A/t^)  model  in  figure  12  and 
(A/t)  model  in  figure  13.  In  both  cases  a  similar 
distribution  in  coda-Q  is  observed  between  the  different 
components.  This  indicates  that  each  component  is  providing 
a  measure  of  similar  earth  properties.  Coda-Q  is  seen  to 
range  between  400  and  2800  for  both  the  (A/t^)  and  (A/t) 
model  results. 

Shown  in  figures  14  and  15  are  the  (A/t^)  and  (A/t) 
model  results  respectively  for  coda-Q  versus  frequency  by 
event.  These  displays  have  a  similar  distribution  and 
variance  in  coda-Q  values  as  observed  in  the  component 
sorted  results.  If  large  differences  were  observed  between 
the  two  shallow  explosions  (events  1  and  2)  and  the  deep 
earthquake,  it  could  imply  near-source  scattering  effects. 
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Figure  12-  (A/t^)  model  coda-Q  versus  frequency  results  by 

common  component  for  all  events  and  stations.  The 
least-square  best  fit  line  is  shown.  Symbols  represent  the 
three  events;  Event  1  (+) ,  Event  2  (x)  and  Event  3  (o) . 
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Figure  13.  (A/t)  model  coda-Q  versus  frequency  results  by 

common  component  for  all  events  and  stations.  The 
least-square  best  fit  line  is  shown.  Symbols  represent  the 
three  events:  Event  l  (+) ,  Event  2  (x)  and  Event  3  (o) . 
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Figure  14.  (A/t^)  model  coda-Q  versus  frequency  results  by 

common  event  for  all  stations  and  components.  The 
least-square  best  fit  line  is  shown.  Symbols  represent  the 
3  components:  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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Figure  15.  (A/t)  model  coda-Q  versus  frequency  results  by 

common  event  for  all  stations  and  components.  The 
least-square  best  fit  line  is  shown.  Symbols  represent  the 
3  components:  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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No  significant  difference  is  observed  suggesting  that 
near-source  scattering  effects,  if  present,  are  poorly 
resolved. 

Coda-Q  versus  frequency  results  plotted  by  common 
receiving  station  are  given  in  figure  16  for  the  (A/t^) 
model  and  figure  17  for  the  (A/t)  model.  For  common 
receiver  sorting  the  coda-Q  solutions  are  observed  to 
cluster  along  linear  trends  indicating  a  strong  dependence 
on  frequency  and  near- receiver  earth  properties.  The  large 
decrease  in  variability  in  both  the  (A/t^)  and  (A/t)  model 
displays  with  the  common  receiver  sorting  indicates  that 
near-receiver  scattering  is  the  dominant  contributor  to 
teleseismic  coda  formation. 

The  (A/t^)  model  results  in  figure  16  for  stations 
RSCP  and  RSNY  between  1  and  2  Hz  appears  to  have  higher 
coda-Q  values  than  the  other  receiving  stations.  The 
variance  is  also  greater  at  RSCP  indicating  a  lower 
signal-to-noise  ratio  for  this  station.  RSCP  and  RSNY  also 
appear  to  have  anomalous  slopes  of  coda-Q  with  frequency 
compared  to  RSON,  RSNT  and  RSSD.  To  test  if  the  slope  of 
the  coda-Q  with  frequency  relationship  at  station  RSCP  is 
similar  to  that  at  RSNY  the  poor  data  at  RSCP  were  edited 
and  the  least-square  analysis  repeated.  The  edited  and 
unedited  data  are  shown  in  figure  18.  The  least-square  best 
fit  line  for  station  RSCP  after  editing  is  equal  to  that  at 
station  RSNY. 

A  similar  argument  for  the  RSCP  and  RSNY  station 
is  more  difficult  for  the  (A/t)  model  results  in  figure  17. 
This  is  due  to  the  very  high  variability  in  these  displays 
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Figure  16.  (A/t^)  model  coda-Q  versus  frequency  results  by 

common  receiving  station  for  all  events  and  components.  The 
least-square  best  fit  line  is  shown.  Symbols  represent  the 
3  components;  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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Figure  17 •  (A/t)  model  coda-Q  versus  frequency  results  by 

common  receiving  station  for  all  events  and  components.  The 
least-square  best  fit  line  is  shown.  Symbols  represent  the 
3  components:  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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Figure  18.  Edit  test  of  (A/t^)  model  coda-Q  solutions  at 
station  RSCP.  Trend  at  station  RSCP  is  similar  to  station 
RSNY  after  editing  poor  data. 
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which  make  the  trends  harder  to  interpret. 

The  observed  range  of  coda-Q  values  for  both  the 
(A/t^)  and  (A/t)  model  results  are  very  similar  in  all  the 
displays.  The  fact  that  both  models  produce  similar  coda-Q 
estimates  is  in  agreement  with  the  synthetic  tests  in  paper 
2.  The  (A/t^)  model  appears  to  produce  results  with  lower 
variance  compared  to  the  (A/t)  model.  The  best  fit  lines  in 
figure  17  for  the  (A/t)  model  results  do  not  follow  the 
data  trend  as  well  as  in  the  (A/t^)  model  results  in  figure 
16.  The  large  amount  of  scatter  in  the  data  points  is 
corrupting  the  least-square  fitting  process. 

The  above  observations  suggest  that  qualitatively  the 
(A/t^)  model  fits  the  data  better  than  the  (A/t)  model,  but 
is  that  supported  with  quantitative  evidence?  One  measure 
of  which  model  fits  the  data  better  is  the  number  of  stable 
inversion  solutions  obtained  by  each  model  for  the  same 
data.  Another  measure  of  which  model  best  fits  the  data  is 
the  sum  of  the  square  of  the  residuals  (referred  to  as  SSR 
in  the  text  that  follows) .  The  SSR  is  a  measure  of  how 
closely  the  data  points  follow  the  model  curve.  Table  3 
gives  the  number  of  stable  solutions  and  SSR  for  both 
models  on  a  representative  suite  of  the  data. 

To  effectively  use  the  number  of  stable  inversions 
and  the  SSR  as  a  measure  of  best  model,  it  must  be  assumed 
that  a  stable  inversion  always  implies  a  proper  fit  and 
therefore  a  good  coda-Q  er'^.imate.  Based  on  the  high 
variances  and  evidence  that  the  (A/t)  model  appears  to  be 
finding  more  solutions  above  5  Hz  where  the  signal-to-noise 
ratio  is  poor  (figures  16  and  17,  station  RSCP)  this 
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assumption  is  probably  not  valid.  It  is  therefore  difficult 
to  draw  any  conclusions  from  table  3.  The  small  magnitude 
of  the  differences  and  the  fact  that  no  consistent  trend  is 
observed  suggests  that  both  models  fit  the  data  equally 
well. 


Table  3. — Quantitative  comparison  of  the  (A/t^)  and  (A/t) 
model  fits  to  the  data,  shown  are  the  number  of  stable 
solutions  and  the  SSR  for  a  representative  subset  of  data 


NUMBER  OF  SOLUTIONS  BY  EVENT  (0-10  Hz,  Z-Component  at  RSON) 
MODEL  EVENT  1  EVENT  2  EVENT  3 


(A/t2) 

28 

29 

36 

(A/t) 

22 

26 

27 

NUMBER  OF 

SOLUTIONS 

BY  COMPONENT  (0-10 

Hz,  Event  2  at  RSNT) 

MODEL 

Z-COMPONENT 

R-COMPONENT  T-COMPONENT 

(A/t2) 

34 

30 

40 

(A/t) 

23 

40 

53 

NUMBER  OF 

SOLUTIONS 

BY  STATION 

(0-10  Hz 

;,  Event  3,  Z-comp.) 

MODEL 

RSON 

RSNT 

RSSD 

"’.SCP  RSNY 

(A/t2) 

36 

6 

23 

20  15 

(A/t) 

27 

11 

26 

19  18 

SUM  OF  THE  SQUARE  OF  THE  RESIDUALS  BY  EVENT  (Station  RSON) 

EVENT  1  EVENT  2  EVENT  3 


FREQUENCY 

(A/t2)  (A/t) 

(A/t2)  (A/t) 

(A/t2)  (A/t) 

.8 

626 

651 

307 

306 

611 

640 

1.4 

861 

834 

1262 

1254 

1032 

1051 

2.1 

1551 

1536 

1707 

1726 

543 

518 

3.0 

624 

624 

361 

338 

561 

561 

Another  means  of  quantifying  which  of  the  models  is  a 
better  fit  to  the  data  that  is  also  sensitive  to  the 
quality  of  the  final  results  is  the  least-square  best  fit 
line.  The  assumption  is  made  that  the  coda-Q  values  follow 
a  linear  trend  increasing  with  frequency.  The  coda-Q 
solutions  from  both  models  appear  to  suggest  this  type  of 
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relationship.  The  standard  deviation  of  the  y-intercept  and 
slope,  the  X-Y  correlation  coefficient  and  residuals  from 
the  least-square  best  fit  line  analysis  is  given  in  table 
4.  The  standard  deviations  of  the  slope  and  y-intercepts 
are  comparable  for  all  stations  between  the  two  models. 
However,  the  X-Y  correlation  coefficients  are  higher  and 
the  residuals  lower  for  the  (A/t^)  model  results,  except  at 
station  RSSD.  These  results  are  in  agreement  with  the 
visual  appearance  of  plotted  coda-Q  values. 


Table  4. — The  least-square  best  fit  line  analysis 


(A/t2)  MODEL 

STANDARD  DEVIATION 

X-Y 

STATION 

Y-INTERCEPT 

SLOPE 

CORRELATION 

RESIDUALS 

RSON 

50 

14 

.81 

346 

RSNT 

80 

24 

.68 

454 

RSSD 

79 

33 

.60 

427 

RSCP 

125 

51 

.23 

536 

RSNY 

125 

69 

.73 

501 

(A/t)  MODEL 

STANDARD  DEVIATION 

X-Y 

STATION 

Y-INTERCEPT 

SLOPE 

CORRELATION 

RESIDUALS 

RSON 

53 

16 

.64 

437 

RSNT 

67 

21 

.61 

491 

RSSD 

55 

23 

.69 

389 

RSCP 

105 

33 

.28 

550 

RSNY 

148 

68 

.20 

756 

The  major  difference  between  the  (A/t^)  and  (A/t) 
model  results  appears  to  be  a  larger  degree  of  variance  in 
the  (A/t)  model  result  displays.  This  suggests  that  the 
(A/t^)  model  fits  the  data  better  and  produced  more 
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reliable  coda-Q  estimates  than  the  (A/t)  model.  A  possible 
explanation  for  this  can  be  seen  in  figure  2.  Because  the 
time  dependent  term  is  not  well  separated  from  the  time 
independent  ^erm  the  inversion  explains  some  of  the  later 
part  of  the  coda  with  the  former.  Note  how  the  model  curve 
drawn  through  the  latter  part  of  the  coda  in  the  lower 
display  in  figure  2  does  not  follow  the  data  as  well  as  in 
the  upper  display.  However  when  the  time  dependent  and 
independent  terms  are  summed  and  the  resulting  model  curves 
displayed  over  the  data  (figure  5)  there  is  good  agreement 
in  both  cases. 

The  inability  of  the  (A/t)  model  to  fit  steep  decays 
in  the  early  part  of  the  coda  suggests  that  the  (A/t)  model 
is  explaining  mainly  the  later  part  of  the  coda.  This 
splitting  of  the  later  coda  energy  between  the  two  terms 
may  result  in  biased  coda-Q  estimates  which  could  explain 
the  high  variances  in  the  displays.  Comparison  of  the 
(A/t^)  and  (A/t)  models  in  paper  2  using  synthetics 
supports  this  interpretation.  Based  on  this  limited 
comparison  the  (A/t^)  model  appears  to  describe  coda 
with  a  large  multiple  energy  component  at  teleseismic 
distances  better  than  the  (A/t)  model. 

Because  the  (A/t^)  model  appears  to  produce  more 
consistent  results  they  will  be  used  in  the  following 
quantitative  comparisons  to  other  studies  of  earth  Q.  To 
quantify  the  frequency  dependence  of  coda-Q,  f*^  is  computed 
from  the  least-square  best  fit  line  for  the  (A/t^)  model 
results  sorted  by  receiving  station.  The  value  n  is 
computed  from  results  at  two  different  frequencies  by  the 
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following  relation: 

(  Q/Qo  )  =  (  f/fo  )" 

STATION  =  RSNT  RSON  RSSD  RSNY  RSCP 

n  =  .56  .52  .54  .72  .16(.73) 

Stations  RSON,  RSNT  and  RSSD  have  very  similar  coda-Q 
versus  frequency  relations.  Stations  RSCP  and  RSNY  appear 
to  have  anomalous  slopes  that  could  be  due  to  the  lower 
signal-to-noise  ratios  observed  at  these  stations.  Based  on 
the  discussion  for  figure  16,  the  results  for  RSCP  may  be 
in  error  due  to  low  signal-to-noise  ratio.  A  value  of  n 
».73  similar  to  RSNY  was  computed  for  station  RSCP  when  the 
poor  data  was  edited  prior  to  least-square  fitting.  Values 
of  slope  reported  by  Jin  et  al.,  (1985)  ranged  between 
.46-. 61  for  local  events  in  old  (stable)  and  young  (active) 
oceanic  areas  respectively.  The  slopes  computed  in  this 
study  are  in  qualitative  agreement  with  theirs. 

Figure  19  is  a  comparison  between  coda-Q  results  of 
this  study  and  those  of  Singh  and  Herrmann  (1975)  from 
local  and  near  regional  coda  analysis.  A  similar 
relationship  is  observed  between  the  coda-Q  results  from 
this  study  and  their  independent  results  using  short 
source-to-receiver  separations.  With  short  source  to 
receiver  separations  scattering  processes  near  the  source 
or  receiver  provide  measures  of  similar  earth  properties. 
Agreement  between  local  and  teleseismic  coda-Q  results 
supports  the  conclusion  that  teleseismic  coda-Q  estimates 
are  controlled  by  near  receiver  scattering  and  provide  a 
measure  of  earth  properties  near  the  receiving  station. 

The  frequency  dependence  of  teleseismic  coda-Q 
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Figure  19.  (A/t^)  model  coda-Q  versus  frequency  solutions 

by  common  receiver.  Plotted  is  the  least-square  best  fit  to 
the  solutions  from  this  study,  event  1  (*),  2  (o)  and  3 
(X) .  The  local  coda-Q  results  from  Singh  and  Herrmann 
(1983)  are  also  shown  (•). 
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estimates  are  compared  to  those  from  local  coda-Q  analysis 
in  figure  20.  Good  agreement  is  observed  between  both  the 
absolute  magnitude  and  slope  of  the  coda-Q  with  frequency 
relations. 

CONCLUSIONS 

The  GLI  approach  described  in  this  paper  is  an 
alternative  method  for  model  fitting  in  coda  analysis.  It 
has  the  advantage  of  using  the  whole  coda  where  signal  is 
sufficiently  above  noise  rather  than  small  subsets  of  the 
coda. 

The  (A/t2+B)exp(-Ct)  model  produced  more  consistent 
results  than  a  (A/t+B)exp(-Ct)  model.  Application  of  the 
combined  single  plus  triple  body  wave  backscattering  model 
(Gao  et  al.,  1983b)  with  the  (A/t^+B) exp (-Ct)  form  provided 
a  means  of  obtaining  stable  estimates  of  coda-Q  when 
multiple  scattering  effects  are  present.  Coda-Q  determined 
from  fitting  the  single  scattering  model  alone  when 
multiple  scattering  effects  are  present  would  produce 
biased  solutions.  A  Monte  Carlo  error  analysis  indicates 
that  reliable  results  are  obtainable  with  the  GLI  coda 
analysis  approach  for  levels  of  random  noise  observed  in 
teleseismic  data.  The  details  of  the  Monte  Carlo  test  and 
the  bias  test  are  the  subject  of  paper  2. 

For  the  teleseismic  P-wave  codas  analyzed  in  this 
study  a  model  which  accounts  for  both  single  and  multiple 
scattering  effects  was  found  to  be  necessary.  Due  to  the 
strong  multiple  scattering  effects  the  observational  data 
could  not  be  properly  fit  with  the  single  scattering  model 
alone.  The  results  of  a  test  to  convert  the  A  and  B  model 
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FREQUENCY  in  Hr 


Figure  20.  The  frequency  dependence  of  coda-Q  for  local 
earthquakes  as  summarized  by  Aki  (1980a, b).  The  lowest 
frequency  data  are  those  from  fundamental  surface  waves. 
The  teleseismic  coda-Q  results  from  this  study  are  plotted 
for  comparison. 


parameters  into  earth  turbidity  and  source  spectral 
estimates  were  encouraging,  a  more  detailed  study  of  this 
is  the  subject  of  paper  3. 

The  coda-Q  values  predicted  in  this  study  are  found 
to  be  in  agreement  with  other  estimates  in  the  literature. 
Coda-Q  values  increased  with  frequency  in  all  cases  in  the 
range  from  approximately  800  to  1600  from  1  to  4  Hz.  Slopes 
of  the  best  fit  lines  ranged  between  f*52  ^^d  f*^^.  The 
slope  at  station  RSCP  computed  from  the  least-square  best 
fit  line  was  questionable  due  to  low  signal-to-noise  ratio, 
but  was  found  to  be  closer  to  that  of  station  RSNY 
when  editing  is  performed  prior  to  least-square  fitting. 
Higher  coda-Q  values  of  1600  at  i  Hz  were  observed  at 
stations  RSCP  and  RSNY,  compared  to  1100  at  stations  RSON, 
RSNT  and  RSSD. 

A  strong  relationship  between  coda-Q  and  near¬ 
receiver  scattering  was  observed.  The  coda-Q  displays 
have  a  large  variance  except  when  the  results  are  plotted 
by  common  receiving  station.  No  near-source  scattering 
differences  were  observed  between  shallow  explosions  and  a 
deep  earthquake.  In  fact,  almost  identical  coda-Q  estimates 
were  obtained  for  these  two  cases.  There  is  good  agreement 
between  the  coda-Q  results  of  this  study  and  other 
published  Q  values  from  local  and  regional  work.  These 
agreements  support  the  conclusion  that  the  coda-Q  values 
estimated  from  teleseismic  coda  using  the  methods  described 
in  this  paper  provide  information  about  earth  properties 
near  the  receiver.  Near-source  scattering  effects  were  not 
found  to  be  resolvable  in  this  study. 
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PAPER  II 


SUMMARY 

Two-terro  models  of  the  form  (A/t^+B) exp(~Ct)  and 
(A/t+B)exp(-Ct)  where  found  to  characterize  both  the  single 
and  multiple  scattering  effects  in  coda.  Both  of  these 
models  will  estimate  coda-Q  values  that  are  equal  within 
error  but  the  A  and  B  model  parameters  are  very  different. 
Inversion  stability  may  be  useful  in  determining  which  of 
these  models  is  the  appropriate  one  to  use  for  different 
data  types.  Instabilities  occurred  because  the 
(A/t+B) exp(-ct)  model  does  not  adequately  describe  steep 
decays  at  early  times.  The  A  and  B  model  solutions  were 
affected  most  by  this  deficiency.  The  two-term  model 
parameter  estimates  were  found  to  be  consistent  with 
respect  to  random  noise,  parameter  magnitude  and  type  of 
multiple  scattered  energy  present  in  the  data. 

Strongly  biased  solutions  were  obtained  when  the 
scattering  model  used  in  the  analysis  differed 
significantly  from  the  type  of  scattering  present  in  the 
data.  Bias  is  defined  as  the  percent  difference  between  a 
known  solution  and  one  computed  from  a  coda  analysis 
inversion.  When  the  single  scattering  model  was  fit  to  data 
with  multiple  scattering,  coda-Q  values  were  biased  by  as 
much  as  280%  and  the  product  of  turbidity  and  source  term 
by  267%.  The  bias  in  coda-Q  was  29%,  turbidity  55%  and 
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source  term  170%,  when  the  (A/t^+B) exp(-Ct)  model  was  fit 
to  data  which  contained  only  single  scattering. 

The  two-term  models  were  fit  to  data  containing 
different  orders  of  multiple  scattering  and  the  bias 
decreased  as  additional  orders  of  multiple  scattering,  up 
to  4th  order,  were  added  to  the  input  data.  The  bias  to 
coda-Q  in  this  case  was  9%,  turbidity  14%  and  source  term 
110%.  This  bias  is  similar  in  magnitude  to  that  introduced 
by  10-20%  random  noise.  Similar  to  the  effects  of  random 
noise,  the  bias  due  to  fitting  inappropriate  models  was 
found  to  affect  the  source  term  parameter  the  most  and 
coda-Q  the  least.  No  significant  frequency  dependent  bias 
was  observed. 

Tests  performed  with  the  energy- flux  model  showed 
that  the  two-teino  models  produced  similar  Q  estimates  when 
intrinsic  Q  effects  dominated  over  scattering  Q  effects 
(scattering  Q  is  large) .  Agreement  between  these  models 
implies  that  two-term  models  can  be  used  in  place  of  the 
single  scattering  model  for  coda-Q  estimation  when  multiple 
scattering  effects  are  present  in  the  data. 

INTRODUCTION 

This  is  the  second  of  three  papers  which  address  the 
issue  of  how  to  deal  with  multiple  scattered  energy  in  coda 
analysis.  The  first  paper  (Soroka  and  Stump,  1991)  deals 
with  coda-Q  estimation  and  will  be  referred  to  as  paper  1 
in  the  following  text.  The  third  paper  (Soroka,  Stump  and 
Dainty,  1991)  referred  to  as  paper  3  in  the  text,  deals 
with  earth  turbidity  and  source  spectral  estimation  in  coda 
analysis.  In  this  paper  a  study  using  synthetic  data  is 
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undertaken  to  determine  if  simple  two-term  mooei-.  to 

used  in  coda  analysis  to  describe  single  plus  mu  in i pie 
scattered  coda  energy.  The  effects  of  different  typer; 
of  scattered  energy,  coda  duration  and  random  noise  on  coda 
analysis  solutions  are  described.  Different  scattering 
theories  are  compared  to  better  understand  their 
similarities,  differences  and  limitations  to  help  justify 
the  application  of  two-term  models.  The  results  of  this 
study  also  demonstrate  the  effect  an  incorrect  model  could 
have  on  coda  analysis  estimates  when  it  is  difficult  to 
determine  the  correct  model  to  use. 

The  analysis  of  coda,  the  random  energy  that 
immediately  follows  a  body  wave  arrival,  has  been  shown  to 
provide  useful  estimates  of  earth  and  source  properties 
such  as  Q,  turbidity  and  source  spectrums  (Aki  and  Chouet, 
1975;  Herraiz  and  Espinosa,  1988;  Rautian  and  Khalturin, 
1978;  paper  1).  A  coda  analysis  involves  fitting  a  model  to 
observational  data  and  then  computing  earth  and  source 
properties  from  the  model  parameters. 

Because  it  is  not  always  possible  to  record  data 
close  to  an  event,  coda  analysis  work  has  begun  to  be 
extended  to  teleseismic  distances  where  codas  are  dominated 
by  multiple  scattering  effects.  The  approach  used  by 
some  authors  has  been  to  fit  the  single  scattering  model  in 
the  same  manner  as  with  local  and  regional  studies,  even 
though  multiple  scattering  effects  were  suspected  j’*  the 
data.  The  approach  used  successfully  in  paper  1  is  to  fit  a 
two-term  model  to  better  account  for  the  multiple  scattered 
energy.  A  two-terra  model  approach  would  be  equally 
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applicable  at  local  and  regional  distances  when  multiple 
scattering  effects  are  present  in  the  data. 

Selection  of  the  appropriate  model  to  fit  to 
observational  data  in  a  coda  analysis  is  a  very  important 
consideration.  In  paper  1  the  (A/t^+B) exp(-Ct)  and 
(A/t+B) exp(-Ct)  models  hereafter  referred  to  as  the  (A/t^) 
and  (A/t)  models  respectively  were  found  to  fit  the  spectra 
of  coda  dominated  by  multiple  scattered  energy.  The  (A/t^) 
model  was  found  to  produce  coda-Q  estimates  with  lower 
variances  than  the  (A/t)  model  in  that  study.  Figure  l  is  a 
representative  fit  of  the  (A/t^)  model  to  time  varying 
teleseismic  P-wave  coda  spectra.  In  the  upper  display  the 
time  dependent  and  time  independent  terms  are  plotted 
separately  and  in  the  lower  display  the  sum  of  the  two 
terms  is  plotted.  There  is  good  agreement  between  the  model 
and  data.  The  time  dependent  term  is  necessary  to  describe 
the  early  single  scattering  effects  and  the  time 
independent  term  is  necessary  to  describe  the  later 
multiple  scattering  effects.  The  (A/t)  model  was  found  to 
produce  a  similar  match  but  was  deficient  in  modeling  steep 
decays  at  early  times. 

In  this  second  paper  an  objective  is  to  quantify  the 
appropriateness  of  using  simplified  (A/t^)  or  (A/t)  models 
in  coda  analysis  so  that  multiple  scattered  energy  can  be 
included  in  the  analysis.  To  better  understand  the 
appropriateness  of  using  these  simple  two-term  models, 
comparisons  are  made  between  them  and  the  commonly  used 
single  scattering  model  originally  proposed  by  Aki  (1969) 
and  expanded  to  3-D  body  waves  by  Gao  et  al.  (1983b)  and 
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Figure  1.  The  (A/t^+B) exp(-Ct)  model  fit  to  time  varying 
spectral  data  (points) •  In  the  upper  case  the  two  terms  are 
plotted  seperatly  and  in  the  lower  case  their  sum  is 
plotted.  The  data  is  for  Event  3  at  station  RSON,  2  Hz. 


the  energy-flux  model  of  Frankel  and  Wennerberg  (1987) . 

The  approach  that  will  be  used  to  better  understand 
the  bias  issues  is  to  utilize  synthetics  where  model 
parameters  determined  by  inversion  can  be  compared  directly 
with  known  solutions.  Synthetic  tests  give  us  insight  into 
errors  expected  when  obseirvational  data  is  analyzed.  Bias 
is  computed  as  the  percent  difference  between  the  analysis 
estimate  and  the  known  value. 

I  known  -  estimate  | 

bias  =  (  -  )  100 

known 

Variance  (standard  deviation)  by  comparison  is  a  measure  of 
the  amount  a  value  can  be  expected  to  vary  around  the  true 
value.  Variance  gives  the  uncertainty  in  a  solution.  The 
average  from  many  experiments  would  yield  the  true  solution 
if  no  bias  were  present. 

TWO-TERM  MODEL  ASSUMPTIONS 

It  is  common  practice  to  represent  complicated 
geophysical  processes  by  simplified  models  in  the  hope  of 
gaining  insight  into  the  fundamentals  of  the  problem.  The 
wide  use  of  Aki's,  one  term,  single  scattering  model  to 
estimate  earth  Q  is  just  one  example  which  demonstrates  the 
validity  of  this  approach.  Extension  of  the  one  term  single 
scattering  model  to  two-terms  so  that  multiple  scattering 
effects  can  be  used  in  a  coda  analysis  is  also  a 
simplification  of  a  complex  problem. 

In  order  to  map  the  model  parameters  to  more  useful 
earth  and  source  property  estimates  the  (A/t^)  and  (A/t) 
models  must  be  equated  to  an  appropriate  scattering  theory. 
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The  (A/t^)  model  is  the  form  of  first  and  third  order  body 
wave  backscattered  energy  (Gao  et  al.,  1983b).  This  theory 
is  based  on  elastic  wave  propagation  with  the  following 
assumptions: 

1)  The  seismic  energy  emanates  from  a  point  source 

2)  The  source  and  receiver  are  located  at  the  same 

position 

.  3)  The  Born  approximation  is  used  and  so  energy  is  not 

conserved  in  the  system  of  equations 

The  effect  of  the  above  assumptions  is  one  of  the 
issues  of  this  study.  These  assumptions  are  minimized  in 
the  teleseismic  problem  if  the  scattered  waves  which  form 
the  coda  can  be  regarded  as  the  superposition  of  secondary 
waves  which  are  the  sum  of  many  small  independent  events 
resulting  from  the  scattering  interactions  (Aki,  1969) .  If 
the  scattering  occurs  close  to  the  receiver  the  effect  of 
the  coincident  source-receiver  and  point  source  assumptions 
are  therefore  minimized.  Strong  evidence  was  found  in  paper 
1  to  suggest  that  near  receiver  scattering  is  responsible 
for  most  coda  generation  at  teleseismic  distances.  The 
effect  of  the  Born  approximation  is  more  difficult  to 
♦  resolve;  however,  coda-Q  estimates  in  paper  1  are  in 

agreement  with  other  independent  Q  estimates.  Comparison  of 
the  (A/t^) ,  (A/t)  and  energy-flux  models  later  in  this 
paper  also  shows  that  these  assumptions  may  not  present  a 
significant  problem,  at  least  to  the  coda-Q  estimates. 

The  above  assumptions  are  not  as  serious  at  short 
source-receiver  separations  as  they  are  at  large 
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separations.  In  fact,  to  minimize  the  effect  of  the  above 
assumptions  most  coda  analysis  work  has  involved  local  and 
regional  distance  data,  with  small  source-receiver 
separations  the  effects  of  both  source-receiver  coincidence 
and  the  point  source  assumptions  are  thus  minimized.  The 
fact  that  reasonable  parameter  estimates  are  obtained  that 
agree  with  other  independent  measures  of  the  same 
parameters  suggests  that  the  Born  approximation  is  also  not 
a  significant  problem. 

The  (A/t)  model  represents  an  alternative  two-term 
model  for  describing  coda  which  contains  multiple  scattered 
energy.  Like  the  (A/t^)  model,  the  (A/t)  model  must  be 
equated  to  an  appropriate  scattering  theory  in  order  to 
interpret  the  model  parameters.  The  (A/t)  model  was  shown 
to  be  the  correct  form  of  scattered  surface  waves  (Gao  et 
al.,  1983a).  Because  of  the  supporting  evidence  that  coda 
generation  is  due  to  body  wave  scattered  energy  (paper  1) , 
surface  waves  have  been  ruled  out  as  a  significant 
contributor  to  coda  generation. 

The  (A/t)  model  also  has  the  correct  form  of  a  plane 
wave  striking  a  layer  of  random  scatterers  plus  P  to  Lg 
(trapped  shear)  waves  (paper  3) .  This  alternative  theory 
provides  a  means  of  mapping  the  model  parameters  of  the 
(A/t)  model  to  earth  and  source  estimates.  This  plane  wave 
theory  is  more  applicable  to  coda  at  teleseismic  than  local 
distances.  Because  time  is  measured  relative  to  the  direct 
arrival  in  this  case  and  there  are  no  assumptions  regarding 
source-receiver  coincidence  this  theory  is  appealing. 
However,  higher  order  body  wave  scattering  is  assumed  to  be 
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negligible  in  this  case. 

In  paper  l  it  was  shown  that  both  the  (A/t^)  and 
(A/t)  models  produced  reasonable  fits  to  coda  spectral  data 
which  contained  multiple  scattered  energy.  The  two  theories 
will  be  compared  quantitatively  in  the  inversion  results 
section  to  determine  the  degree  of  bias  that  would  result 
if  the  wrong  model  was  applied  in  the  coda  analysis. 

MODEL  OVERVIEW 

Because  the  (A/t^)  model  appeared  to  produce  better 
coda~Q  results  in  paper  1,  it  will  be  fit  to  the  different 
data  sets  throughout  this  study  and  will  serve  as  a  basis 
for  comparison.  In  order  to  determine  the  effect  of  the 
different  tests  on  the  inversion  solutions  the  (A/t^)  model 
parameters  will  be  converted  to  earth  and  source  property 
estimates  using  the  Gao  et  al.  (1983a  and  b)  theory.  Single 
scattering  is  defined  as  one  scattering  interaction  between 
the  source  and  receiver.  Double  scattering  occurs  when  a 
wave  undergoes  a  scattering  interaction  and  the  resulting 
scattered  wave  undergoes  a  second  scattering  interaction 
before  reaching  the  receiver.  Third  and  fourth  order 
scattering  are  defined  in  similar  fashion.  Expressions  for 
the  first  four  orders  of  scattering  are: 


SINGLE 

2  T  S 

DOUBLE 

Pl(f,t)  = 

V  t2 

2.46  t2 

exp(-27rft/Q) 

s 

(1) 

P2(f/t)  = 

t 

—  exp(-27rft/Q) 

(2) 

TRIPLE 

P3(f.t)  = 

.716  T-^ 

V  S  exp (-27rft/Q) 

(3) 
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QUADRUPLE 


.51  v2  s  t  exp(-27rft/Q)  (4) 


P4(f.t)  = 

Where:  P(f,t)  =  Spectral  power  with  frequency  and  time 

T  =  Turbidity 

S  =  Source  term 

V  =  Average  medium  velocity 

Q  =  Coda-Q 

Coda-Q  is  the  inverse  of  apparent  attenuation  and 
includes  the  effects  of  both  intrinsic  attenuation  and 
attenuation  due  to  scattering.  Turbidity  is  a  measure  of  a 
medium's  ability  to  initiate  scattering  and  the  source  term 
is  an  estimate  of  the  source  power  spectrum  (adapted  from 
Gao  et  al .  1983a  and  b) .  Use  of  the  single  scattering 
model  alone  does  not  allow  for  the  separation  of  turbidity 
and  source  term.  However,  with  the  (A/t^)  model,  which 
corresponds  to  single  plus  triple  order  scattering  in  the 
abc'^e  equations,  separation  of  the  source  and  turbidity 
terms  is  possible. 

The  four  types  of  scattering,  equations  1-4,  are 
displayed  in  figure  2  for  comparison.  Single  scattering  is 
seen  to  dominate  the  early  part  of  the  coda  and  decays 
quickly  with  time.  Double  scattering  begins  with  lower 
magnitude  than  the  single  scattering  model  but  decays  more 
slowly.  After  a  short  time  the  double  model  curve  rises 
above  the  single  model  curve  indicating  the  dominance  of 
double  over  single  scattering  in  the  coda  decay  curve.  The 
triple  scattering  curve  is  a  straight  line  on  this  log 
display  because  there  is  no  time  dependence  other  than  the 
exponential  term  in  this  model.  At  long  times  third  order 
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Figure  2.  Comparison  of  different  types  of  backscattered 
body  wave  models.  Shown  are  the  single  through  guadruple 
models  and  the  combined  single  plus  triple  (A/t^)  model. 
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scattering  rises  above  the  single  and  double  mode]  curves. 
Fourth  order  scattering  has  very  low  magnitude  at  early 
time  but  increases  in  amplitude  until  it  rises  above  all 
the  other  models  at  longer  times. 

In  summary  the  early  part  of  the  coda  decay  curve 
will  be  dominated  by  single  scattering  interactions.  For  a 
short  transition  time  double  scattering  will  dominate 
followed  by  triple  and  possibly  quadruple  order  scattering 
very  late  in  the  coda.  The  single  plus  triple  scattering 
model  is  also  plotted  in  figure  2  for  comparison.  It  is 
difficult  to  see  because  it  is  on  or  near  the  single 
scattering  curve  at  early  time,  then  the  double  scattering 
curve  for  a  short  time  and  finally  the  triple  scattering 
curve  at  later  time.  The  single  plus  triple  model  appears 
to  represent  a  good  average  model  of  all  these  higher  order 
scattering  effects. 

Because  the  (A/t^)  model,  equation  5,  which  includes 
first  and  third  order  scattering  fits  the  data  well  a  model 
that  also  included  second  order  scattering,  equation  6, 
might  be  even  better. 

A 

(  “  +  B  )  exp(-C  t)  (5) 

t2 

A  D 

(  -  +  -  +  B  )  exp(-C  t)  (6) 

t2  t 

Tests  using  a  single  plus  double  plus  triple  scattering 
model  produced  very  unstable  results.  The  extra  free 
parameter  "D”  in  equation  (6)  was  not  resolvable  due  to  the 
short  transition  time  in  which  this  2nd  order  scattering 
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term  Is  the  dominant  contributor. 

Since  simple  two-term  models  may  not  represent  all 
types  of  scattering  present  in  the  data  it  is  important  to 
understand  how  the  contributions  of  these  other  forms  of 
scattered  energy  will  affect  the  inversion  solutions.  The 
bias  that  results  from  ignoring  double  and  quadruple 
scattering  will  be  addressed  latter  in  this  paper. 

PARAMETER  SENSITIVITY  SYNTHETICS 
By  combining  the  four  models  discussed  above 
(equations  1-4)  five  synthetics  were  created  for  input  to 
the  coda  analysis  program.  The  five  models  are  listed  below 
as  well  as  the  abbreviation  that  will  be  used  throughout 
the  remainder  of  this  paper. 

SYNTHETIC  DATA  TYPES  ABBREVIATION 


1) 

Single 

(single) 

2) 

Single  + 

Double 

(1+2) 

3) 

Single  + 

Double  + 

Triple 

(1+2+3) 

4) 

Single  + 

Double  + 

Triple  +  Quadruple 

(1+2+3+4) 

5) 

Single  + 

Triple 

(1+3) 

For  coda  analysis  it  is  the  combined  effects  of  the 
individual  models  displayed  in  figure  2  that  is  important. 
Analysis  of  these  synthetics  can  provide  insight  on  the 
resolution  of  these  different  orders  of  scattering  by  the 
inversion.  All  five  models  include  single  scattering 
effects  but  differ  in  the  type  of  higher  order  scattering 
added.  The  five  models  are  compared  in  figure  3.  The  single 
scattering  model  is  the  lowest  curve  at  long  times  on  the 
display.  This  is  due  to  the  rapid  decay  of  energy  with  time 
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Figure  3.  Comparison  of  the  synthetic  model  data  analyzed 
in  this  study.  All  the  curves  contain  single  scattering  but 
different  orders  of  multiple  scattering. 


when  only  first  order  scattering  is  present.  As  additional 
orders  of  scattering  are  added  the  decay  with  time  becomes 
increasingly  slower.  The  (A/t^)  model  is  equivalent  to  the 
1+3  model  and  can  be  seen  to  fall  between  the  1+2  and  1+2+3 
models.  All  the  models  converge  toward  the  single 
scattering  model  at  early  times. 

INVERSION  PROCEDURES 

The  five  synthetic  data  sets  described  above  were 
used  as  input  to  a  coda  analysis  using  Generalized  Linear 
Inversion  (GLI)  techniques.  Approximately  1%  random  noise 
was  added  to  each  data  set  to  stabilize  the  inversion. 

Both  the  single  scattering  model,  equation  (1),  of  Gao  et 
al.  (1983b)  and  the  (A/t^)  model,  equations  (1)  plus  (3), 
were  fit  to  the  synthetic  data.  For  the  single  model  fit 
coda-Q  (Q)  and  the  product  of  turbidity  and  source  term 
(TS)  are  computed.  With  the  (A/t^)  model  fit  it  was 
possible  to  separate  the  turbidity  and  source  term 
parameters.  I  .  this  case  turbidity  (T) ,  source  term  (S)  and 
coda-Q  (Q)  estimates  are  reported. 

It  was  observed  that  the  inversions  would  converge 
toward  a  reasonable  solution  for  most  input  models. 

That  is,  the  sum  of  the  square  of  the  residuals  (referred 
to  as  SSR  in  the  text  that  follows)  between  model  and  data 
would  decrease  in  a  systematic  manner  toward  a  minima.  Low 
SSR's  meant  that  the  model  and  data  were  in  good  agreement, 
high  SSR's  meant  poor  agreement.  If  the  input  model  was 
significantly  different  from  the  model  being  fit  low  SSR's 
may  not  be  possible.  In  such  cases  unrealistic  parameters 
(negative  values)  could  be  computed  with  eventual  failure 


of  the  inversion  if  the  criteria  to  end  the  inversion  was 
not  satisfied.  For  example,  when  fitting  the  (A/t^) 
two-term  model  to  single  scattering  input  the  inversion 
would  be  able  to  reduce  the  SSR's  up  to  a  point.  To  achieve 
a  better  fit  the  amplitude  of  the  time  independent  term 
would  have  to  go  to  zero  so  that  only  the  time  dependent 
term  would  be  fit  to  the  input.  This  suggests  that  the 
inversion  approach,  is  model  sensitive  and  a  stable 
solution  is  unlikely  if  the  observational  data  are 
significantly  different  from  the  model  being  fit.  Rather 
than  reporting  that  the  inversion  failed  or  that  a 
parameter  value  of  zero  was  obtained,  the  model  fit  with 
the  lowest  SSR  before  the  instability  occurred  is  reported. 

INVERSION  RESULTS 
THE  (A/t+B)exp(-Ct)  MODEL  TESTS 

In  paper  1  both  the  (A/t^)  and  (A/t)  models  were 
found  to  describe  coda  dominated  with  multiple  scattered 
energy.  Although  the  (A/t^)  model  solutions  appeared  better 
because  of  lower  variances  in  the  displays  there  was  still 
room  for  debate.  Because  of  the  difficulty  in  determining 
the  appropriate  model  to  use  in  a  coda  analysis  these 
models  are  compared  here  to  better  understand  their 
similarities  and  differences. 

The  (A/t)  model  is  the  correct  form  of  scattered 
surface  waves  and  scattered  energy  from  a  plane  wave 
striking  a  near  surface  layer.  The  time  dependent  term  is 
the  correct  form  for  single  scattered  events  which  result 
when  the  plane  wave  strikes  random  scatterers.  The  time 
independent  term  represents  P  to  Lg  energy  trapped  in  the 
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layer  under  the  receiver.  The  model  parameters  A  and  B  in 
(A/t)  model  are  interpreted  differently  depending  upon 
whether  the  surface  or  plane  wave  model  is  assumed.  The  C 
term  from  which  coda-Q  is  computed  is  however  the  same  for 
all  the  models,  including  the  (A/t^)  model. 

As  a  way  of  comparing  the  two  models  synthetic  data 
are  generated  with  the  (A/t)  model  and  the  (A/t^)  model  is 
fit  to  that  data.  The  A  and  B  terms  from  the  input  model 
are  compared  to  the  inversion  solutions  to  determine  what 
effect  the  different  time  dependence  of  the  single 
scattering  term  has  on  the  model  solutions.  The  known 
coda-Q  value  is  also  cc;xipared  with  the  coda-Q  inversion 
solution  to  determine  the  bias  that  would  result  from  using 
these  different  models. 

Table  1  gives  the  (A/t^)  model  coda-Q  results  for  a 
range  of  different  input  model  values.  In  each  case  the 
correct  value  of  coda-Q  was  determined,  within  error,  by 
fitting  the  (A/t^)  model.  Similar  results  were  obtained 
when  different  values  of  A  and  B  were  used  in  the  input 
model . 


Table  1. — The  1+3  model  Coda-Q  solutions  for  a  range  of 

input  model  Q’s 


KNOWN  CODA-Q 

INVERSION  CODA-Q 

%  DIFFERENCE 

1000 

998 

<  1% 

800 

799 

<  1% 

600 

599 

<  1% 

400 

400 

0% 

200 

200 

0% 

Shown  in  table  2  are  the  results  of  a  test  to 
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determine  the  effect  of  limiting  the  number  of  data  points 
in  the  GLI  analysis.  The  coda*-Q  solutions  appear  to  be 
relatively  robust  in  all  cases  tested.  In  the  extreme  case 
when  only  10  samples  were  used  only  12%  difference  was 
observed  between  the  known  and  inversion  coda-Q  values. 
This  could  be  attributed  to  the  1%  random  noise  added  to 
the  synthetic  data  prior  to  analysis.  These  tests  suggest 
that  very  similar  coda-Q  solutions  would  be  obtained  by 
fitting  either  the  (A/t)  or  (A/t^)  models  to  observational 
data.  In  paper  1  these  two  models  were  fit  to  a  suite  of 
teleseismic  coda  and  the  range  of  coda-Q  values  obtained 
were  the  same  in  both  cases.  The  results  from  fitting  the 
(A/t2)  model  had  less  variance  than  the  (A/t)  model 
suggesting  the  (A/t^)  model  fit  the  characteristics  of  the 
data  better.  No  evidence  was  found  in  the  tests  on 
synthetic  to  suggest  a  preference  of  one  model  over  the 
other  in  coda-Q  estimation. 


Table  2. — (A/t2)  model  coda-Q  solutions  as  a  function  of 
number  of  analysis  points,  input  data  generated  with  the 
(A/t)  model,  2  Hz  case  with  1%  random  noise 


Number  of 
samples 

KNOWN 

CODA-Q 

INVERSION 

CODA-Q 

PERCENT 

DIFFERENCE 

160 

1000 

999 

<  1% 

100 

1000 

998 

<  1% 

75 

1000 

979 

2% 

40 

1000 

977 

2% 

20 

1000 

910 

9% 

10 

1000 

1123 

12% 

The  A  and  B  model  parameters  will  now  be  examined  to 
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complete  the  comparison  of  the  (A/t^)  and  (A/t)  models.  The 
input  model  A  and  B  parameters  are  compared  to  those 
determined  by  fitting  the  (A/t^)  model  in  table  3.  The 
results  show  significant  differences  between  the  A  and  B 
values.  The  low  SSR's  indicate  that  a  good  fit  of  the 
(A/t^)  model  to  the  data  was  obtained  in  each  case.  These 
good  fits  were  obtained  by  decreasing  the  amplitude  of  the 
time  dependent  (A)  term  in  the  model  and  forcing  the  B  term 
and  exponential  to  explain  the  data.  In  other  words,  the 
(A/t^)  term  is  inadequate  to  describe  (A/t)  term  effects 
early  in  the  coda  data  but  the  later  coda  is  described  well 
with  the  B  term.  Higher  levels  of  random  noise  or 
decreasing  the  number  of  analysis  points  did  not  change 
this  situation. 

Table  3. — Comparison  of  the  input  model  and  inversion 
solution  A  and  B  model  parameters 


INPUT 

MODEL 

FINAL  SOLUTION 

SSR 

A 

B 

A 

B 

5.0 

0.0875 

1.3E-10 

5.1 

8.3 

10.0 

0.35 

1.9E-17 

10.5 

8.3 

50.0 

8.75 

7.9E-11 

59.3 

8.3 

100.0 

35.0 

1.2E-10 

136.2 

8.3 

The  results  of  this  test  suggest  that  low  residuals 
can  be  obtained  by  fitting  either  of  these  two  models  to 
observational  data.  Dramatically  different  values  of  the  A 
and  B  model  parameters  will  be  obtained  depending  upon 
which  model  is  used.  Coda-Q  estimates  will  be  very  similar 
regardless  of  which  model  is  used.  Selection  of  the 
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appropriate  model  to  use  is  however  important  if  the  A  and 
B  model  parameters  are  to  be  properly  converted  to  earth 
and  source  property  estimates.  The  subject  of  converting 
the  A  and  B  model  parameters  into  earth  and  source 
properties  is  discussed  in  paper  3. 

(A/t2+B)exp(-Ct)  MODEL  SENSITIVITY  TESTS 

The  characteristics  of  data  used  in  coda  analysis  can 
be  highly  variable.  For  the  results  of  a  coda  analysis  to 
be  useful  they  must  be  reasonably  consistent  with  changes 
in  signal-to-noise  ratio,  duration  of  the  coda  and  types  of 
scattered  energy  present.  The  following  tests  are  designed 
to  quantify  the  effect  these  types  of  changes  will  have  on 
parameter  estimates  when  fitting  the  (A/t^)  model  to  coda 
spectra . 

The  following  results  were  obtained  by  fituing  the 
(A/t^)  model  to  the  5  synthetic  data  sets  described  above. 
The  sensitivity  of  the  inversions  to  the  order  of 
scattering  present  in  the  input  data  is  quantified.  Table  4 
gives  the  coda-Q  (Q) ,  turbidity  (T)  and  source  term  (S) 
inversion  solutions  for  the  different  input  models.  The 
corresponding  bias  or  percent  difference  between  the  known 
answer  and  the  inversion  result  and  the  SSR  are  also  given. 
An  SSR  value  not  equal  to  zero  for  the  1+3  model  solution 
is  due  to  the  random  noise  added  to  the  synthetics.  This 
random  noise  effect  will  be  similar  for  each  of  the  models 
tested.  Based  on  the  SSR  a  similar  quality  fit  was  made  to 
the  1+2  and  1+2+3  synthetic  data.  The  large  increase  in  SSR 
when  the  single  and  1+2+3+4  synthetics  are  used  suggests 
that  the  characteristics  of  these  data  are  not  well 
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represented  by  the  (A/t^)  model. 


Table  4. — Inversion  solutions  and  bias  from  fitting  the 
(A/t2)  model  to  data  with  different  types  of  multiple 
scattering  present,  2  Hz  case 


MODEL  KNOWN  MODEL  USED  TO  GENERATE  INPUT 


PARAMETER 

SOLUTION 

1 

1+2 

1+3 

1+2+3 

1+2+3+4 

S 

(bias) 

100000 

282547 

(183%) 

263853 

(164%) 

100781 

(1%) 

174500 

(75%) 

90952 

(9%) 

T 

(bias) 

.01 

.004 

(57%) 

.007 

(26%) 

.01 

(0%) 

.01 

(0%) 

.014 

(40%) 

Q 

(bias) 

1000 

777 

(22%) 

764 

(24%) 

1000 

(0%) 

798 

(20%) 

1210 

(21%) 

SSR 

513 

51 

13 

69 

212 

In  the  next  series  of  tests  the  sensitivity  of  the 
coda  analysis  solutions  to  the  magnitude  of  each  model 
parameter  is  studied  over  their  commonly  observed  range. 
The  values  of  Q,  T  and  S  cover  the  ranges  observed  in 
papers  1  and  3.  The  results  in  table  5  show  the  effect  of 
varying  the  number  of  data  points  used  in  the  analysis. 


Table  5.— Number  of  analysis  points  test.  The  (A/t^)  model 
is  fit  to  the  1+2+3  synthetic  data  at  2  Hz 


MODEL 

NUMBER  OF 

DATA  POINTS  USED 

PARAMETER 

known 

160 

100 

50 

30 

10 

S 

100000 

174499 

146118 

104151 

81844 

385405 

(bias) 

(74%) 

(46%) 

(4%) 

(18%) 

(285%) 

T 

.01 

.01 

.011 

.014 

.017 

.002 

(bias) 

(0%) 

(10%) 

(40%) 

(70%) 

(78%) 

Q 

(bias) 

100 

98 

95 

90 

87 

201 

(2%) 

(5%) 

(10%) 

(13%) 

(101%) 

SSR 

69 

40 

13 

28 

32 
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The  biases  to  S,  Q  and  T  were  all  observed  to  decrease  with 
increasing  number  of  analysis  points.  The  source  term  had 
the  most  bias  with  the  least  number  of  analysis  points.  The 
bias  however  was  a  minimum  with  50  points  rather  than  with 
the  maximum  points  as  with  Q  and  turbidity. 

The  effect  of  Q  magnitude  on  solution  bias  is  shown 
in  table  6.  Varying  Q  has  no  effect  on  the  source  term  and 
turbidity  estimates.  Solution  bias  in  coda-Q  is  observed  to 
increase  with  increasing  Q  values.  The  observed  range  of 
coda-Q  in  paper  1  was  between  400  and  2000.  Values  100  and 
5000  were  included  to  test  the  extreme  cases. 

A  similar  test  but  with  the  1+2  synthetic  as  input 
produced  similar  results  to  those  shown  in  table  6:  no 
change  in  T  or  S  bias  and  Q  bias  decreasing  with  decreasing 
Q.  These  results  suggest  that  a  coda  analysis  which  fits 
the  (A/t^)  model  to  data  is  more  robust  against  bias  when 
coda-Q  values  are  low. 

Table  6. — Q  magnitude  test,  inversion  results  and  bias 
from  fitting  the  (A/t^)  model,  2  Hz,  160  points  case 

MODEL  Q  VALUE  USED  IN  1+2+3  INPUT  MODEL 


PARAMETER 

known 

100 

1000 

3000 

5000 

S 

100000 

174500 

174500 

174500 

174500 

(bias) 

(75%) 

(75%) 

(75%) 

(75%) 

T 

.01 

.01 

.01 

.01 

.01 

(bias) 

(0%) 

(0%) 

(0%) 

(0%) 

Q 

— — 

98 

798 

1707 

2210 

■m 

(bias) 

(2%) 

(20%) 

(43%) 

(56%) 

SSR 

69 

69 

69 

69 

Changes 

in  bias  due  to 

varying  the 

magnitude  of 

turbidity  are  shown  in  table  7.  Turbidity  controls  the 
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crossover  point  from  single  to  multiple  scattering  effects 
in  the  (A/t^)  model.  Lower  values  of  turbidity  result  in 
greater  separation  of  single  and  multiple  scattering 
effects  by  delaying  the  crossover  point  in  time.  Values  of 
turbidity  were  observed  to  range  from  approximately  .01  to 
.1  in  paper  3  but  values  less  than  .01  are  published  in  the 
literature.  Turbidity  magnitude  influences  the  S,  Q  and  T 
estimates.  Less  bias  appears  to  occur  with  lower  values  of 
T,  although  there  appear  to  be  exceptions  to  this  general 
observation. 


TABLE  7. — Turbidity  magnitude  test,  shown  is  the  (A/t^) 
model  solutions  and  bias,  2  Hz,  160  data  point  case 


MODEL 

PARAMETER 

T  VALUE  USED  IN  THE  1+2+3  INPUT  MODEL 
known  .0001  .001  .01  .05 

S 

100000 

125209 

117465 

174500 

351453 

(bias) 

(25%) 

(17%) 

(75%) 

(251%) 

T 

... 

.00008 

.0009 

.01 

.035 

(bias) 

(20%) 

(10%) 

(0%) 

(30%) 

Q 

500 

513 

662 

444 

484 

(bias) 

(3%) 

(32%) 

(11%) 

(3%) 

SSR 

13 

92 

69 

33 

Changes  in  bias  due  to  varying  the  magnitude  of  the 
source  term  are  shown  in  table  8 .  Source  term  magnitude 
does  not  influence  the  bias  to  the  S,  Q  or  T  estimates 
significantly. 

A  test  of  the  frequency  dependence  of  the  bias  is 
shown  in  table  9.  The  (A/t^)  model  is  fit  to  synthetic  data 
generated  with  the  1+2  model  at  1,  3  and  10  Hz.  There  is 
essentially  no  change  in  bias  to  the  source  term  and 
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turbidity  results.  The  coda-Q  estimates  improve  with 
increasing  frequency.  Since  frequency  is  in  the  numerator 
of  the  exponential  term  and  coda-Q  :  -  in  the  denominator, 
an  inverse  relationship  should  exist  between  these  two 
parameters.  Based  on  the  previous  testing  of  Q  magnitude, 
bias  to  Q  should  decrease  with  increasing  frequency  and 
there  should  be  no  effect  to  the  T  and  S  estimates. 


Table  8. — Source  term  magnitude  test,  shown  is  the  (A/t^) 
model  solutions  and  bias,  2  Hz,  160  point  case 


MODEL  S  VALUE  USED  IN  THE  1+2+3  INPUT  MODEL 


PARAMETER 

known 

1000 

10000 

100000 

500000 

S 

1827 

17450 

174500 

872497 

(bias) 

(82%) 

(75%) 

(75%) 

(74%) 

T 

.01 

.0096 

.01 

.01 

.01 

(bias) 

(4%) 

(0%) 

(0%) 

(0%) 

Q 

500 

478 

444 

444 

444 

(bias) 

(4%) 

(11%) 

(11%) 

(11%) 

SSR 

13 

92 

69 

33 

Table  9. — Frequency  test,  the  (A/t^)  model  solutions  at 
1,  3  and  10  Hz  with  the  1+2  model  used  as  input 


MODEL 

KNOWN 

FREQUENCY 

PARAMETER 

SOLUTION 

IHZ 

3HZ 

lOHz 

S 

100000 

263853 

263853 

263853 

(bias) 

(164%) 

(164%) 

(164%) 

T 

.01 

.0075 

.0075 

.0075 

(bias) 

(25%) 

(25%) 

(25%) 

Q 

1000 

618 

829 

942 

(bias) 

(38%) 

(17%) 

(6%) 

SSR 

51 

51 

51 
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To  sum  up,  parameter  sensitivity  tests  have  shown 
that  bias  varies  depending  upon  the  magnitude  of  Q  and  T 
but  not  S.  The  observed  bias  is  less  when  Q  and  T  are  lower 
in  magnitude  (100,  .001)  and  greater  when  Q  and  T  are 
larger  (1000,  .05).  Decreasing  Q  (or  increasing  frequency) 
produces  faster  coda  decay  rates  which  allow  for  better 
model  parameter  estimates.  Small  values  of  T  cause  the 
single  to  multiple  scattering  transition  to  occur  at  a 
later  time;  the  additional  single  scattering  information 
improves  model  parameter  resolution.  Bias  in  the  estimates 
was  also  found  to  depend  on  the  number  of  data  points  or 
duration  of  the  coda.  The  more  data  available  to  fit 
the  model,  the  less  bias  there  is  in  the  inversion 
estimates . 

MONTE  CARLO  ERROR  ANALYSIS 

A  Monte  Carlo  error  analysis  was  performed  to 
determine  the  effects  of  random  noise  on  the  (A/t^)  model 
solutions.  Both  the  variance  and  bias  in  the  presence  of 
different  levels  of  random  noise  were  quantified.  The  test 
consisted  of  100  coda  analysis  runs  with  the  same  input 
data  but  with  different  sets  of  random  noise  of  a  constant 
level  added  to  the  input  data  each  time.  A  statistical 
analysis  on  the  Q,  T  and  S  solutions  from  the  lOO  runs  were 
used  to  calculate  the  bias  and  variance  in  each  of  the 
model  parameters. 

Two  tests  on  the  effects  of  random  noise  were  made. 

In  the  first  test  the  1+3  model  was  used  to  generate  the 
input  and  the  (A/t^)  model  is  fit  to  the  data.  In  this  case 
the  model  to  be  fit  and  the  model  used  to  generate  the 
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synthetic  data  are  the  same.  Figure  4  gives  the  expected 
variance  (percent  error  or  standard  deviation)  in  coda-Q, 
turbidity  and  source  term  as  a  function  of  random  noise. 

Also  shown  is  a  range  of  random  noise  observed  in  paper  1. 

Based  on  the  sensitivity  tests  and  the  fact  that  high 
magnitude  values  of  Q  and  T  were  used  in  this  first  test, 
the  results  in  figure  4  represents  a  worst  case. 

In  the  second  test  (figure  5)  the  1+2+3  model  was  « 

used  to  generate  input  with  low  magnitude  values  of  Q 
and  T,  As  expected  from  the  parameter  sensitivity  tests, 
the  variances  in  the  model  parameter  estimates  were 
significantly  reduced  with  lower  coda-Q  and  turbidity 
values.  The  stronger  decay  of  coda  energy  with  time  and 
better  separation  of  single  and  multiple  coda  information 
allow  for  more  robust  inversions.  Note  that  a  significant 
scale  change  from  that  used  in  figure  4  was  required  to 
make  figure  5. 

Table  10  gives  bias  as  a  function  of  random  noise  for 
the  coda-Q,  turbidity  and  source  term  for  both  cases.  In 
the  first  case  the  input  model  is  the  same  as  the  model 
being  fit  and  values  of  Q  (1000)  and  T  (.1)  are  at  the  high 
end  of  the  typical  range.  Bias  increases  with  increasing 
random  noise  magnitude.  Q  and  T  are  relatively  unaffected 
by  even  20%  random  noise;  however,  S  is  strongly  affected. 

In  the  second  case  the  input  model  was  different  than  the 
model  being  fit  and  values  of  Q  (500)  and  T  (.01)  are 
closer  to  the  low  end  of  their  range.  The  bias  in  this  case 
is  due  primarily  to  the  difference  in  scattering 
characteristics  between  the  input  and  model  being  fit  and 
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MONTE  CARLO  ERROR  ANALYSIS 


OBSERVED 
DATA  NOISE  RANGE 


SOURCE  TERM 


TURBIDITY 


CODA-Q 


%  RANDOM  NOISE 


Figure  4.  Monte  Carlo  error  analysis  results  with  high 
magnitude  parameter  values.  Shown  is  the  variance  (standard 
deviation)  as  percent  error  from  the  true  value  with 
different  levels  of  random  noise  in  the  data.  The  range  of 
random  noise  observed  in  teleseismic  data  is  also  shown. 


MONTE  CARLO  ERROR  ANALYSIS 


» 


Figure  5.  Monte  Carlo  error  analysis  results  with  low 
magnitude  parameter  values.  Shown  is  the  variance  (standard 
deviation)  as  percent  error  from  the  true  value  with 
different  levels  of  random  noise  in  the  data. 
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not  due  to  the  random  noise.  Addition  of  random  noise 
actually  improves  the  Q  estimates  slightly  in  this  test. 
Bias  to  T  and  S  remains  relatively  constant  except  with  the 
40%  random  noise  case  which  causes  the  bias  to  increase. 


Table  10. — Bias  computed  as  percent  difference  between 
the  known  answer  and  the  Monte  Carlo  test  mean 


MODEL 

PARAMETER 

1 

RANDOM 

10 

NOISE 

20 

40% 

1+3  INPUT  MODEL 

S 

0% 

30% 

200% 

2000% 

T 

0% 

3% 

10% 

115% 

Q 

0% 

1% 

3% 

10% 

1+2+3  INPUT  MODEL 

S 

52% 

53% 

56% 

73% 

T 

16% 

16% 

17% 

22% 

Q 

27% 

27% 

26% 

21% 

The  parameter  sensitivity  tests  and  Konte  Carlo 
analysis  show  that  fitting  the  (A/t^)  model  to  coda  with 
varying  random  noise,  duration  and  orders  of  scattering 
produces  reasonably  robust  solutions.  The  bias  problem  is 
less  with  smaller  values  of  Q  and  turbidity.  Table  11 
compares  the  Monte  Carlo  analysis  results  to  the  results 
from  fitting  the  (A/t^)  model  to  synthetic  data  with  up  to 
4th  order  scattering.  Only  the  results  for  low  values  of  Q 
and  turbidity  and  10  and  20%  random  noise  are  shown.  The 
magnitude  of  the  bias  due  to  using  the  wrong  model  is 
equivalent  to  that  due  to  10%  to  20%  random  noise.  The 
results  in  table  11  should  be  compared  with  those  in  table 
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4,  which  were  determined  with  large  Q  and  turbidity  values. 


Table  11. — Comparison  of  the  bias  due  to  random  noise  and 
the  effect  of  fitting  an  inappropriate  model,  2  Hz  case 


MODEL  MODEL  USED  TO  GENERATE  INPUT  RANDOM  NOISE 


PARAMETER 

1 

1+2 

1+2+3 

1+2+3+4 

10% 

20% 

S  1000 

2% 

4% 

3% 

18% 

30% 

200% 

T  .001 

2% 

4% 

4% 

16% 

3% 

10% 

Q  200 

1% 

5% 

5% 

8% 

1% 

3% 

SINGLE  SCATTERING  MODEL  TEST 

The  bias  due  to  fitting  the  single  scattering  model 
(equation  1)  to  data  which  include  multiple  scattering 
information  is  included  in  this  study  because  other  authors 
have  demonstrated  its  inability  to  correctly  model  data 
with  multiple  scattering  effects  (Frankel  and 
Wennerberg,  1987) .  The  effect  of  incorrectly  fitting  the 
(A/t^)  model  to  data  composed  of  only  single  scattering 
information  is  also  of  interest.  Table  12  gives  the 
inversion  results  obtained  by  fitting  the  single  scattering 
model  to  the  4  synthetic  data  sets  which  contain  up  to  3rd 
order  scattering.  The  inversion  with  the  single  model 
failed  when  4th  order  scattering  was  included. 

Because  there  is  only  one  amplitude  term  involved  in 
fitting  the  single  model  alone,  it  is  not  possible  to 
separate  the  turbidity  and  source  terms.  The  product  of 
turbidity  and  source  term  is  therefore  reported  in  table 
12.  The  bias  is  consistently  higher  than  when  the  (A/t^) 
model  was  fit  to  these  same  data  (table  4)  .  This  is  due  to 
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the  fact  that  higher  order  scattering  effects  present  in 
the  data  are  not  being  properly  accounted  for  when  the 
single  model  is  fit  to  these  data.  These  results  imply  that 
coda-Q  solutions  computed  by  fitting  the  single  scattering 
model  can  be  strongly  biased  if  the  coda  data  being  fit  are 
dominated  by  multiple  scattered  energy. 


Table  12. — Inversion  results  when  a  single  scattering 
model  is  fit  to  the  synthetic  input  data,  2  Hz  case 


MODEL 

KNOWN 

MODEL  1 

USED  TO 

GENERATE 

INPUT 

PARAMETER 

SOLUTION 

SINGLE 

1+2 

1+2+3 

1+3 

TS 

1000 

1005 

2907 

3667 

1965 

(bias) 

(1%) 

(191%) 

(267%) 

(97%) 

Q 

1000 

1000 

1618 

2904 

3804 

(bias) 

(0%) 

(62%) 

(190%) 

(280%) 

ENERGY-FLUX  MODEL  TEST 


The  Energy-Flux  model  has  been  used  by  Frankel  and 
Wennerberg  (1987)  and  Langston  (1989)  to  model  coda.  The 
Energy-Flux  model  is  based  on  the  diffusion  equation  which 
accounts  for  multiple  scattering  and  provides  a  means  to 
separate  scattering  Q  and  intrinsic  Q  effects.  While  the 
(A/t^)  model  allows  for  the  separation  of  the  source, 
turbidity  and  coda-Q  terms  it  cannot  differentiate 
scattering  Q  from  intrinsic  Q.  The  Energy-Flux  model  from 
Langston  (1989)  is: 


A(f,t)= 


1  1 

27rftd(-  +  -  ) 

(2Idtd)l/2  Qs  Qi 

-  exp( - 

t  2 


(7) 


-27rft  -27rft 

)  (l-exp(  — ) )  exp( - ) 

2Qs  2Qi 
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Where:  A(f,t)  =  coda  spectral  amplitude 


t(j  =  travel  time  of  P-wave  through  layer 
Qg  =  scattering  Q  parameter 
Qi  =  intrinsic  Q  parameter 

I(j  =  integral  of  direct  wave  squared  velocity 

f  =  frequency 

t  =  time 

As  a  means  of  comparing  the  (A/t^)  model  to  the 
Energy-Flux  models  test  data  were  produced  with  the 
Energy-Flux  model  and  then  analyzed.  The  quality  factor  is 
computed  from  Qi  and  Qs  using  equation  8  for  comparison  to 
the  coda-Q  estimates: 

1  11 

=  -  +  -  (8) 

Quality  factor  Qg 

Table  13  compares  the  Energy-Flux  model  parameters  to 
the  (A/t^)  model  estimates.  At  Q  values  of  1000  the  (A/t^) 
model  coda-Q  was  biased  by  only  2%  from  the  correct  answer 
suggesting  that  the  two  models  are  similar  in  this  case. 
With  Q  values  of  500  a  35%  bias  was  observed  and  at  Q 
values  of  200  a  bias  of  69%  was  observed.  Figure  6  shows 
the  strong  similarity  between  the  (A/t^)  model  fit  to  the 
Energy-Flux  model  data  for  both  a  low  and  high  Q  test. 

Because  the  coda-Q  estimates  in  table  13  appeared  to 
closely  follow  the  values  of  Qi,  two  tests  were  performed 
in  which  Qi  and  Qs  were  varied  separately.  In  the  first 
test  shown  in  table  14,  Qs  is  kept  constant,  at  1000  and  Qi 
is  varied  from  800  to  100.  For  all  values  of  Q  the  coda-Q 
estimates  appear  to  agree  very  well  with  the  quality  factor 
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Figure  6.  Comparison  of  the  energy-flux  model  to  the 
(A/t^+B) exp (-Ct)  model.  The  upper  display  is  with  coda-Q 
equal  to  500  and  the  lower  display  is  with  coda-Q  equal  to 
333.  The  (A/t^)  model  (line)  is  fit  to  data  generated  with 
the  energy  flux  model  plus  0.5%  random  noise  (points). 


computed  from  equation  8. 


Table  13. — Energy-Flux  model  test  with  intrinsic  and 
scattering  Q  equal,  the  (A/t^)  model  results  with  the 
the  Energy-Flux  model  used  as  input,  quality  factor  is 

defined  in  equation  8 


Energy-Flux  INPUT  PARAMETERS 

QS 

200 

500 

1000 

Qi 

200 

500 

1000 

Quality  factor 

100 

250 

500 

GLI  RESULTS 

coda-Q 

169 

337 

532 

(bias) 

(69%) 

(35%) 

(6%) 

Turbidity 

.008 

.016 

.021 

Source  term 

343988 

56512 

18324 

Table  14. — Energy-Flux  model  test  with  intrinsic  and 
scattering  Q  different,  the  1+3  model  inversion 
results  when  the  Energy-Flux  model  is  used  as  input, 
cpaality  factor  is  computed  from  equation  8 


INPUT  PARAMETERS 


Qs 

1000 

1000 

1000 

1000 

Qi 

800 

500 

300 

100 

Quality  factor 

444 

333 

230 

91 

GLI  RESULTS 

coda-Q 

470 

344 

233 

92 

(bias) 

(6%) 

(3%) 

(1%) 

(1%) 

In  the  second  case,  shown  in  table  15,  Qi  is  kept 
constant  at  1000  and  Qs  is  varied  from  800  to  100.  The 
coda-Q  estimates  in  this  case  remain  relatively  constant  at 
approximately  500.  Changing  the  magnitude  of  Qs  does  not 
appear  to  significantly  alter  the  rate  of  decay  in  time  of 
the  Energy-Flux  model.  Although  these  tests  are  limited, 
they  suggest  that  it  is  intrinsic  Q  that  controls  the  decay 
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of  energy  with  time  with  scattering  Q  affectino  the  decay 
to  a  much  lesser  degree  in  the  Energy-Flux  model.  The 
(A/t^)  model  appears  to  agree  closely  with  the  energy-flux 
model  except  in  the  case  were  scattering  Q  is  small 
relative  to  intrinsic  Q. 


Table  15. — Energy-Flux  model  test  with  intrinsic  and 
scattering  Q  different,  the  1+3  model  inversion  results 
when  the  Energy-Flux  model  is  used  as  input,  quality 
factor  is  computed  from  equation  8 


INPUT  PARAMETERS 


Qs 

800 

500 

300 

100 

Qi 

1000 

1000 

1000 

1000 

Quality  factor 

444 

333 

233 

91 

GLI  RESULTS 

coda-Q 

508 

487 

495 

602 

(bias) 

(14%) 

(46%) 

(112%) 

(562% 

CONCLUSIONS 

Two-term  models  were  found  to  describe  coda  spectral 
decay  with  time  and  therefore  appear  useful  in  modeling 
multiple  scattered  energy  in  coda.  The  two  term  models 
consisted  of  a  time  dependent  term  for  the  steep  decays 
early  in  the  coda  and  the  time  independent  term  for  the 
gentle  decays  later  in  the  coda.  The  two  models  studied  in 
this  paper  were:  (A/t^+s) exp(-Ct)  and  (A/t+B) exp(-Ct) . 
Comparison  of  the  two  models  showed  that  similar  coda-Q 
values  would  be  computed  with  either  model.  Large 
differences  were  however  observed  between  the  A  and  B  model 
parameters  and  the  (A/t)  model  appeared  deficient  at 
modeling  steep  spectral  decays.  These  differences  suggest 
that  model  stability  during  inversion  could  help  determine 
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the  most  appropriate  model  to  use  in  a  coda  analysis. 

Parameter  sensitivity  tests  with  the  (A/t^)  model 
showed  that  solution  bias  is  dependent  upon  the  magnitude 
of  Q  and  T  but  not  S.  The  bias  was  less  when  low  Q  (100) 
and  T  (.001)  values  were  used  and  higher  with  large  Q 
(1000)  and  T  (.05)  values.  Bias  in  the  estimates  was  also 
found  to  depend  upon  the  number  of  data  points  or  duration 
of  the  coda.  The  more  data  available  the  less  bias  there 
was  in  the  inversion  estimates.  This  observation  suggestr 
that  better  coda  analysis  estimates  are  possible  if 
two-term  models  which  use  more  of  the  coda  are  used  instead 
of  just  a  single  scattering  model  applied  to  only  the  very 
early  part  of  the  coda. 

Biased  solutions  were  obtained  when  the  models  to  be 
fit  differed  significantly  from  the  scattering  present  in 
the  data.  The  amount  of  bias  was  found  to  depend  on  the 
degree  of  difference.  When  a  single  scattering  model  was 
fit  to  data  with  multiple  scattering  effects  coda-Q  values 
were  biased  by  as  much  as  280%  and  the  product  of  turbidity 
and  source  term  by  as  much  as  267%.  When  the  (A/t^)  model 
was  fit  to  data  consisting  of  only  single  scattering, 
coda-Q  values  are  found  to  be  biased  by  22%,  turbidity  by 
57%  and  source  term  by  183%. 

When  the  (A/t^)  model,  which  includes  multiple 
scattered  effects,  was  fit  to  data  with  different  orders  of 
multiple  scattered  energy  the  amount  of  bias  to  the 
estimates  was  reduced.  The  (A/t^)  model  was  found  to  be  a 
good  average  model  for  data  which  contained  up  to  4th  order 
multiple  scattering.  For  large  values  of  Q  (1000),  T  (.01) 
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and  S  (100000).  the  maximum  bias  in  coda-Q  was  24%, 
turbidity  40%  and  source  term  164%.  By  comparison,  the 
bias  due  to  approximately  20%  random  noise  is  3%  for 
coda-Q,  10%  for  turbidity  and  200%  for  the  source  term. 

When  small  values  of  Q  (200),  T  (.001)  and  S  (1000)  were 
used  the  maximum  biases  were  8%  for  Q,  16%  for  T  and  18% 
for  S. 

The  solution  biases  due  to  fitting  inappropriate 
models  was  found  to  be  greatest  for  the  source  term  and 
least  for  coda-Q.  Solution  bias  due  to  random  noise  was 
also  greatest  for  S  and  least  for  coda-Q.  Coda-Q  is 
therefore  the  most  reliable  estimate  from  a  coda  analysis. 
No  frequency  dependent  bias  was  observed. 

A  marked  reduction  in  stability  of  rhe  inversions  was 
noticed  when  inappropriate  models  were  fit  to  the  data.  The 
analysis  was  found  to  be  model  sensitive.  This  observation 
suggests  that  the  GLI  coda  analysis  technique  is  subject  to 
poor  convergence  or  failure  if  an  unreasonable  model  is 
used. 

When  the  Energy-Flux  model  was  used  to  generate  input 
data  for  the  GLI  coda  analysis  the  solutions  indicated  that 
the  two  models  were  similar  when  scattering  Q  values  were 
high  (1000) .  When  intrinsic  Q  was  high  and  scattering  Q  low 
(100)  the  coda-Q  solutions  differed  from  the  quality  factor 
by  562%  indicating  a  divergence  between  the  two  models.  In 
each  case  coda-Q  was  computed  to  be  larger  than  the  true 
value,  probably  due  to  the  fact  that  higher  order 
scattering  effects  are  included  in  the  Energy-Flux  model. 
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PAPER  III 


SUMMARY 

Two-term  models  of  the  form  (A/t2+B)exp(-Ct)  and  , 

(A/t+B) exp(-Ct)  are  fit  to  time  varying  spectra  of  coda  and 
earth  turbidity  and  source  spectra  estimated.  The  analysis  » 

is  performed  on  teleseismic  P-wave  coda  with  lengths  up  to 
700  seconds.  The  long  durations  suggest  that  the  coda  may 
contain  significant  multiple  scattered  energy.  Separate 
estimates  of  turbidity  and  source  spectrum  were  only 
possible  with  the  (A/t^+B) exp (-Ct)  model  which  assumes  that 
coda  generation  is  due  to  single  and  multiple  scattered 
body  waves.  Turbidity  values  ranged  from  .007  to  .2,  in 
agreement  with  other  published  values.  Turbidity  was  also 
found  to  be  constant  with  freguency.  The  estimated  source 
spectra  had  well  developed  corner  frequencies  and  decay  of 
energy  at  high  frequencies.  Comparison  of  turbidity  and 
source  spectra  products  from  both  models  after 
normalization  showed  them  to  be  equal  within  error.  This 
similarity  suggests  that  the  methods  used  in  this  study  may 
not  be  as  model  sensitive  as  originally  believed.  For  both 
models,  the  fall  off  of  high  frequency  energy  was  observed 
to  be  greater  for  the  explosion  data  than  for  the  deep 
earthquake,  in  agreement  with  other  published  work.  While 
the  results  of  this  feasibility  study  are  encouraging, 
additional  work  on  larger  data  sets  is  necessary  to  better 
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evaluate  the  significance  of  the  turbidity  and  source 
spectra  estimates. 


INTRODUCTION 

This  paper  is  the  last  of  three  which  document 
research  conducted  to  better  understand  how  to  treat 
multiple  scattered  energy  in  a  coda  analysis,  with  primary 
emphasis  on  teleseismic  observations.  Coda  is  the  random 
energy  that  immediately  follows  a  direct  body  wave  arrival. 
In  the  first  paper  (Soroka  and  Stump,  1991) ,  hereafter 
referred  to  as  paper  1,  coda-Q  estimation  using  two-term 
models  and  an  inversion  approach  was  developed.  In  the 
second  paper  (Soroka,  1991),  referred  to  as  paper  2,  the 
results  of  an  error  analysis  study  designed  to  better 
understand  the  applicability  of  simple  two-term  models  in 
coda  analysis  are  given.  Based  on  the  encouraging  results 
in  papers  1  and  2,  the  study  was  expanded  to  determine  if 
additional  earth  and  source  information  could  be  estimated 
from  the  two-term  model  fits  to  the  coda. 

This  third  paper  deals  with  estimating  the  source 
spectrum  and  earth  turbidity  from  the  two-term  model 
parameters  determined  in  paper  1.  Source  spectra  play  an 
important  part  in  the  earthquake-explosion  discrimination 
problem  and  the  coda  method  has  potential  for  producing  a 
more  stable  spectrum  than  by  using  the  direct  arrival 
alone.  Turbidity  is  a  difficult  property  to  estimate  but 
provides  useful  information  about  how  homogeneous  it  is  in 
the  subsurface.  Turbidity  values  have  the.  potential  for 
being  an  indicator  of  tectonic  activity  in  an  area  (Herraiz 
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and  Espinosa,  1987) . 

Because  it  is  not  always  possible  to  record  data 
close  to  an  event,  coda  analysis  work  has  begun  to  be 
extended  to  teleseismic  distances  where  multiple  scattering 
and  crustal  structure  effects  may  be  more  dominant 
contributors  to  the  coda.  The  approach  used  by  some 
researchers  has  been  to  fit  the  single  scattering  model  in 
the  same  manner  as  with  local  and  regional  studies 
(Novelo-Casanova  and  Butler,  1986) .  More  recently  the 
energy  flux  model  has  been  used  to  describe  coda  (Langston, 
1989;  Frankel  and  Wennerberg,  1987).  The  alternative 
approach  developed  in  paper  1  is  to  fit  a  two  term  model  to 
better  account  for  multiple  scattering  and  crustal 
structure  effects  that  occur  later  in  the  coda.  A  two-term 
model  approach  would  also  be  applicable  at  local  and 
regional  distances  when  crustal  structure  and  multiple 
scattering  effects  are  present  in  the  coda.  Use  of 
additional  coda  information  in  a  coda  analysis  could  in 
some  cases  produce  a  better,  more  stable  solution. 

In  current  coda  analysis  a  model  is  fit  to  coda 
spectral  power  as  a  function  of  time.  Then  estimates  of 
earth  and  source  properties  are  computed  from  the  model 
parameters.  Two  important  assumptions  in  this  approach  are 
first  that  a  proper  model  is  being  fit  and  second  that  the 
coda  contains  earth  and  source  information.  Papers  1  and  2 
discuss  the  issue  of  model  choice  in  detail.  There  is 
strong  evidence  in  the  literature  that  coda-Q  estimates  are 
a  measure  of  earth  Q  (Aki  and  Chouet,  1975  and  Herraiz  and 
Espinosa,  1989) .  Useful  information  about  the  source  has 
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also  been  obtained  by  analyzing  coda  (Rautian  and 
Khalturin,  1978) .  In  addition,  Aki  (1982)  found  that  coda 
for  a  particular  event  at  different  receiving  stations  has 
similar  amplitude  and  spectral  characteristics.  These 
observations  suggest  that  the  coda  does  contain  information 
about  the  source  that  initiated  the  seismic  event  and  about 
the  earth  material  through  which  the  seismic  energy  has 
traveled . 

The  principal  earth  information  contained  in  the  coda 
is  Q  and  turbidity.  Earth  Q  is  the  inverse  of  attenuation 
and  is  a  measure  of  the  loss  of  primary  wave  energy  due  to 
scattering  and  absorption.  Local  and  regional  coda  analysis 
work  has  been  aimed  primarily  at  Coda-Q  estimation  (Aki  and 
Chouet,  1975;  Herraiz  and  Espinosa,  1989).  Turbidity  is  a 
measure  of  the  earth's  ability  to  scatter  seismic  energy 
(Dainty  et  al.,  1987).  Work  on  turbidity  is  not  as  abundant 
as  coda-Q  probably  because  of  the  greater  difficulty  in 
making  this  estimate. 

The  principal  source  information  in  the  coda  is  the 
spectrum  of  the  source  wavelet  and  the  source  magnitude. 
Estimates  of  the  source  spectrum  from  coda  have  been  used 
to  calculate  seismic  moments,  corner  frequencies  and  stress 
drops  (Rautian  and  Khalturin,  1978) .  Source  spectra  are 
also  used  in  earthquake-explosion  discrimination  and  event 
magnitude  estimation  studies. 

A  strong  motivation  for  using  the  coda  to  obtain  the 
source  spectrum  is  the  potential  for  producing  a  more 
stable  spectrum  than  if  the  direct  arrival  is  used  alone. 
The  difficulty  lies  in  extracting  this  information  from  the 
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coda  in  a  meaningful  manner.  In  traditional  coda  analysis  a 
one  term,  single  scattering  model  is  fit  to  a  subset  of  the 
coda  shortly  after  the  direct  arrival.  In  this  case  it  is 
not  possible  to  easily  separate  the  source  spectrum  and 
earth  turbidity  information  because  they  appear  as  a 

product.  The  two-term  coda  analysis  approach  described  in  . 

paper  1  models  the  whole  coda  and  has  the  potential  for 

separating  source  spectrum  and  earth  turbidity  information.  ♦ 

This  implies  that  source  spectral  estimation  from  coda  may 
be  possible  when  multiple  scattering  effects  are  included 
in  the  analysis. 

In  this  paper,  estimates  of  source  spectrum  and  earth 
turbidity,  or  their  product  are  determined  from  model 
parameters  after  fitting  the  models  to  a  set  of  teleseismic 
data.  These  results  are  evaluated  to  determine  if  they 
represent  reasonable  estimates  using  general  appearance  and 
comparison  to  other  independent  measures  of  the  estimated 
properties . 


ASSUMPTIONS 

Coda  generation  is  a  complex  process  which  can 
involve  back-  and  forward-scattered  waves,  mode-converted 
waves,  single-  and  multiple- scattered  »/aves  as  well  as 
numerous  other  sources.  Arriving  at  a  theoretical 
description  which  explains  all  these  processes  and  yet  is 
relatively  simple  to  apply  is  a  difficult  if  not  impossible 
problem.  It  is  common  practice  to  represent  complicated 
geophysical  processes  by  simplified  models  in  the  hope  of 
gaining  insight  into  the  fundamentals  of  the  problem.  The 
wide  use  of  Aki's  one-term  model  to  estimate  earth  Q  at 
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local  and  regional  distances  (Herralz  and  Espinosa,  1987) 

Is  one  example  vhlch  demonstrates  the  validity  of  this 
approach.  The  one-term  model  assum.es  an  Infinite  random 
scattering  medium  and  accounts  for  only  single 
backscattered  body  waves. 

Extension  of  the  one-term  single  scattering  model  to 
two  terms  Is  yet  another  simplification  of  the  complex  coda 
problem.  However,  with  the  two-term  model  the  later  part  of 
the  coda  dominated  by  multiple  scattering  effects  can  be 
used  In  the  coda  analysis.  In  this  study  two-term  models  of 
the  form: 

(A/t2+B)exp(-Ct)  and  (A/t+B) exp(-Ct) , 

hereafter  referred  to  as  the  (A/t^)  and  (A/t)  models 
respectively,  are  fit  to  coda  spectra  with  time  (papers  1 
and  2) .  In  both  cases  the  two-term  models  Involve  the  sum 
of  a  time  dependent  and  time  Independent  term  times  an 
exponential  term.  Coda-Q  Is  computed  from  the  exponential 
term  which  Is  the  same  for  both  models.  In  papers  1  and  2 
It  was  shown  that  coda-Q  estimation  Is  not  greatly  affected 
by  which  of  the  two-term  models  Is  fit  to  the  coda  spectra. 

The  difficulty  In  using  simple  models  to  represent 
complex  processes  has  to  do  with  Interpreting  the  model 
parameters  In  terms  of  physical  properties  of  the  process 
under  study.  A  simple  model  which  cannot  be  equated  with  a 
reasonable  physical  process  to  provide  a  way  to  map  the 
model  parameters  to  physical  properties  Is  of  very  little 
use.  On  the  other  hand  If  several  models  are  found  to  fit 
the  daca  or  If  a  ‘•Imple  model  describes  more  than  one 
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physical  process,  some  way  must  be  found  to  choose  among 
the  different  models  and  theory. 

The  (A/t^)  and  (A/t)  models  represent  different 
scattering  processes.  Depending  upon  the  data  set  under 
analysis  one  model  may  be  more  applicable  than  the  other. 

In  paper  2  the  (A/t 2)  model  was  shown  to  represent 
reasonably  single  and  multiple  backscattered  body  waves  in 
the  theory  of  Gao  et  al.  (1983b).  Because  of  a  coincident 
source  and  receiver  assumption  this  model  is  more 
applicable  to  coda  with  small  source-to-receiver 
separations.  At  larger  separations  the  physical 
justification  of  this  model  is  more  difficult  even  if  coda 
generation  is  due  to  scattering  processes  close  to  the 
receiver.  In  this  case  it  is  assumed  that  a  new  point 
source  equal  in  strength  and  spectral  properties  to  the 
incident  wave  is  placed  near  the  receiving  station.  The 
fact  that  reasonable  coda-Q  values  in  close  agreement  to 
the  (A/t)  model  results  were  obtained  in  paper  1  suggests 
that  the  infinite  random  scattering  medium  or  other 
assumptions  imposed  by  the  theory  are  not  violated. 

Ti.'e  (A/t)  model  represents  scattered  body  waves  plus 
P-to-Lg  (trapped  shear)  waves  that  result  when  a  plane  wave 
hits  a  layer  of  random  scatterers  beneath  the  receiver 
(Appendix) .  Because  of  the  plane  wave  assumption  this  model 
appears  more  applicable  to  coda  with  very  large  source-to- 
receiver  separations.  The  (A/t)  model  can  however  have 
applications  at  small  source-to-receiver  separations  if 
suriace  wave  scattering  is  responsible  for  the  coda 
generation  (papers  1  and  2) . 
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Because  it  is  not  clear  how  sensitive  the  desired 
earth  and  source  property  estimates  are  to  the  choice  of 
model  this  work  should  be  looked  on  as  a  feasibility  study. 
The  main  objective  is  to  determine  if  the  coda  analysis 
technique  can  be  extended  to  extract  more  than  just  coda-Q 
from  data  when  two-term  models  are  used.  The  encouraging 
results  of  this  work  should  be  considered  a  beginning  with 
.  more  detailed  studies  performed  on  larger  data  sets  to 

establish  the  stability  of  the  estimates. 

SCATTERING  MODEL  THEORY 

In  paper  1  the  (A/t^)  and  (A/t)  models  were 
successfully  fit  in  a  least-square  manner  to  up  to  500 
seconds  of  teleseismic  coda.  The  difference  between  the  two 
models,  shown  in  equations  1  and  2  below,  is  the  time 
dependence  of  the  first  term.  A,  B  and  C  are  model 
parameters  that  depend  on  frequency.  P(f,t)  is  spectral 
power  of  the  coda  as  a  function  of  both  frequency  and  time. 
Time  begins  with  the  direct  arrival. 

P(f,t)  =  (  A/t2  +  B  )  exp(-C  t)  (1) 

P(f,t)  =  (  A/t  +  B  )  exp(-C  t)  (2) 

Figure  1  gives  the  A,  B  and  C  model  parameters 
determined  by  both  models  for  the  same  event  and  station. 

In  both  cases  the  model  parameters  appear  to  be  measuring 
systematic  changes  in  the  data  that  are  potentially  related 
to  earth  and  source  properties.  The  (A/t^)  model  results 
have  similar  variances  to  the  (A/t)  model  suggesting  that 
both  fit  the  data  well. 

After  the  models  are  fit  to  the  data  and  the  A,  B  and 
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(A/t+B)exp(-Ct)  MODEL 


(A/t2+B)exp(-Ct)  MODEL 


Figure  i.  Comparison  between  the  A,  B  and  C  model 
parameters  determined  by  fitting  the  (A/t^)  and  (A/t) 

models.  Results  are  for  all  components  at  station  RSON  for 
event  l. 
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C  model  parameters  have  been  determined  the  next  step  is  to 
equate  these  parameters  to  earth  and  source  properties  with 
an  appropriate  scattering  theory.  The  3-D  body  wave 
backscattering  model  of  Gao  et  al.  (1983b) ,  which  describes 
coda  spectral  power  variation  with  time  and  frequency  has 
been  used  to  interpret  the  (A/t^)  model  parameters  in 
equation  1.  In  paper  2  the  (A/t^)  model  was  found  to  be  a 
good  average  model  for  single  plus  multiple  backscattered 
body  waves.  The  A/t^  time-dependent  term  has  the  form 
appropriate  for  single  backscattering  and  the 
time- independent  B  term  has  the  form  of  triple  order 
backscattering  as  shown  below; 

ST  ,  -2  TT  f  t 

P(f,t)  =  ( - - - +  .716  S  t3  V  )  exp( - )  (3) 

t2  Q 

where;  S  =  Source  power  spectrum 
T  =  Turbidity 

V  =  Average  earth  velocity 
Q  =  Coda-Q 
f  =  Frequency 
t  =  Time 

The  (A/t)  model,  equation  2,  has  a  form  corresponding 
to  two  different  types  of  scattered  energy.  The  first  is 
scattered  surface  waves  (Gao  et  al.  1983a)  with  a 
coincident  source-receiver  assumption.  Pg  waves  are 
generally  not  considered  a  dominant  contributor  to  coda 
generation  but  Lg  waves  are  in  local  and  regional 
seismograms  (Herrmann,  1980;  Dainty  and  Toksoz,  1990) .  The 
second  is  scattered  energy  that  results  from  a  plane  wave 
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striking  a  layer  of  random  scatterers  (Appendix) .  For  the 
plane  wave  case,  the  time  dependent  term  represents 
scattered  body  wave  energy  of  all  orders  (but  with  single 
scattering  dominant  for  small  values  of  turbidity)  and  the 
time  independent  term  represents  P-to-Lg  (trapped  shear 
wave)  energy  near  the  receiver.  The  plane  wave  model  can  be 
expressed  mathematically  as: 


STh  STV  -2  7rft 

P(f,t)  =  ( - + - )  exp( - )  (4) 

2  t  2  Q 


where  S,  T,  Q  and  V  are  as  defined  above  and  h  is  the 
thickness  of  the  layer. 

For  both  models  the  exponential  term,  from  which 
coda-Q  is  determined,  is  the  same.  In  paper  2  it  was  shown 
that  similar  coda-Q  estimates  would  be  determined  with 
either  of  these  models.  In  paper  1,  where  both  these  models 
were  fit  to  the  same  data,  the  (A/t^)  model  was  found  to 
produce  results  with  better  clustering  and  lower 
variability.  Figure  2  gives  the  coda-Q  results  obtained  by 
fitting  both  the  (A/t)  and  (A/t^)  models  to  the  same  data. 
The  faster  falloff  with  time  of  the  (A/t^)  model  seemed  to 
fit  the  early  part  of  the  coda  better,  suggesting  that  the 
(A/t^)  model  fits  the  data  better. 

In  paper  2  it  was  shown  that  the  A  and  B  model 
parameters  can  be  significantly  different  depending  upon 
which  model  is  fit  to  the  data.  Both  models  are  used  in 
this  study  to  help  establish  the  applicability  of  these 
different  theories  in  coda  analysis  where  multiple 
scattered  energy  and  structure  effects  are  important. 
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Figure  2.  Coda-Q  versus  frequency  for  station  RSNT. 
Comparison  between  the  (A/t^)  and  (A/t)  model  solutions. 
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ANALYSIS  PROCEDURES 


The  A  and  B  model  parameters  determined  with  the 
inversion  procedures  described  in  papers  1  and  2  were  used 
to  generate  the  results  in  this  paper.  A  generalized  linear 
inversion  (GLI)  approach  was  used  to  fit  the  two-term 
models  to  the  data  in  a  least-square  manner.  The  A,  B  and  C 
parameters  are  then  converted  to  earth  and  source  estimates 
using  the  appropriate  theories. 

For  the  (A/t^)  model,  the  3-D  body  wave  theory  of  Gao 
et  al.  (1983b)  is  used.  Equating  equation  1  and  3  gives: 

2  S  T  ,  2  TT  f 

A  =  - ,  B  =  .716  S  t3  V  ,  C  =  - 

V  Q 

Dividing  B  by  A  allows  turbidity  (T)  to  be  computed.  A 
value  of  V=3.5  km/s  is  used  in  the  calculations.  Once  T  is 
determined  S  is  computed  from  A.  The  same  value  for  S  can 
also  be  calculated  from  B. 

For  the  (A/t)  model,  the  plane  wave  theory  (appendix) 
is  used.  Equating  equations  2  and  4  gives: 

STh  STV  2  Tt  f. 

2  '  2  '  Q 

To  compute  the  product  of  T  and  S  in  this  case  A  and  B  are 
summed  and  ST  calculated  by  assuming  values  for  V  and  h 
(V=3.5  km/s  and  h=35  km).  When  ST  is  computed  using  only  A 
or  B,  the  solutions  are  similar  in  shape  but  considerably 
different  in  magnitude.  It  is  not  clear  that  turbidity 
should  be  the  same  for  the  A  and  B  terms  since  different 
types  of  scattering  are  involved. 
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As  the  above  parametrization  indicates  the  (A/t^) 
model  leads  to  separation  of  the  source  spectrum  (S)  and 
earth  turbidity  (T) .  In  the  (A/t)  model,  separation  of 
source  spectra  and  earth  turbidity  is  not  possible  unless 
an  independent  measure  of  either  the  source  spectrum  or 
turbidity  is  available. 

The  teleseismic  P-wave  data  analyzed  in  this  study 
consisted  of  two  near-surface  explosions  and  one  deep 
earthquake  to  contrast  near-source  scattering  effects.  The 
shallow  explosions  may  contain  near-source  crustal 
scattering  effects  but  the  earthquake  does  not.  Each  of  the 
three  events  were  digitally  recorded  at  five  3-component 
North  American  Regional  Seismic  Test  Network  (RSTN) 
stations.  All  five  receiving  stations  were  used  to  contrast 
near-receiver  scattering  effects.  Table  1  lists  the 
characteristics  of  the  data  set  used  in  this  study.  For 
more  detail  the  reader  is  referred  to  paper  1. 

Two  inversions  were  performed  on  the  three  events 
recorded  at  the  five  3 -component  North  American  RSTN 
stations.  In  one  case  the  (A/t^)  model  was  fit  to  the  data; 
in  the  other  the  (A/t)  model  was  fit.  Estimates  of  coda-Q 
from  the  two  inversions  are  reported  in  paper  1.  The  A  and 
B  model  parameters  computed  in  both  inversions  were 
converted  to  earth  and  source  estimates  us in  the 
appropriate  equations  described  above.  These  solutions  are 
compared  and  discussed  in  detail  below. 

RESULTS 

Three  types  of  composite  displays  have  been  assembled 
from  the  results.  First  all  solutions  are  sorted  by  common 
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component  to  determine  if  similar  or  different  information 
has  been  recorded  in  the  vertical,  radial  and  transverse 
directions.  Second  the  solutions  are  sorted  by  common  event 
to  determine  the  significance  of  near-source  scattering  in 
coda  generation.  Third  the  solutions  are  sorted  by  common 
receiving  station  to  emphasize  near-receiver  effects. 


Table  1. — Event  and  station  information 


EVENT  LOCATIONS 

EVENT 

TYPE 

LOCATION 

1 

EXPLOSION 

EASTERN  KAZAKH,  U.S.S.R. 

2 

EXPLOSION 

EASTERN  KAZAKH,  U.S.S.R. 

3 

EARTHQUAKE 

N.  ARGENTINA-CHILE 

STATION  LOCATIONS 

STATION 

LOCATION 

RSCP 

CUMBERLAND  PLATEAU,  U.S.A. 

RSNY 

NEW  YORK,  U.S.A. 

RSSD 

SOUTH  DAKOTA,  U.S.A. 

RSON 

ONTARIO,  CANADA 

RSNT 

NORTHWEST  TERRITORY,  CANADA 

SOURCE-TO  RECEIVER  DISTANCES  IN  DEGREES 

STATION 


EVENT 

MAGNITUDE 

DEPTH 

RSNT 

RSON 

RSNY 

RSSD 

RSCP 

1 

6.1 

NEAR  SURFACE 

68 

79 

83 

86 

94 

2 

6.1 

NEAR  SURFACE 

68 

79 

83 

86 

94 

3 

5.8 

559  Km 

98 

82 

72 

80 

66 

(A/t2+B)exp(-Ct)  MODEL  RESULTS 

In  this  case  it  is  possible  to  separate  the  source 
term  from  turbidity.  These  two  estimates  will  be  presented 
in  turn  beginning  with  turbidity,  which  is  a  measuie  of  the 
medium's  ability  to  scatter.  An  average  crustal  velocity  of 
3.5  km/s  was  used  in  the  calculations. 

The  turbidity  versus  frequency  results  are  plotted  by 
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common  component  in  Figure  3.  A  similar  distribution  is 
observed  in  each  case  with  the  data  values  falling 
approximately  in  the  range  from  .007  to  .2.  Even  though  the 
scatter  is  high  in  these  displays  the  similar  range  in  the 
results  for  each  component  suggests  they  are  measuring 
similar  earth  properties. 

Turbidity  versus  frequency  results  are  plotted  by 
common  event  in  Figure  4.  These  displays  again  show  a 
similar  distribution  of  the  turbidity  values  in  the  range 
from  .007  to  0.2.  Turbidity  estimates  do  not  appear 
significantly  differenc  between  events. 

Turbidity  versus  frequency  results  are  plotted  by 
receiving  station  in  Figure  5.  The  reduction  of  scatter  in 
these  displays  is  not  as  apparent  as  with  the  coda-Q 
results  (paper  1) ,  possibly  because  of  the  higher  variances 
predicted  for  the  turbidity  solutions  by  the  Monte  Carlo 
error  analysis  (paper  2) .  Note  that  only  the  receiving 
station  results  show  significant  differences  in  the  level 
of  the  turbidity,  consistent  with  the  idea  of  near-receiver 
scattering.  The  turbidity  solutions  are  relatively  constant 
over  the  frequency  range  analyzed.  These  displays  suggest 
that  scattering  will  occur  equally  at  all  frequencies 
analyzed  in  this  study.  This  could  be  interpreted  as  an 
indication  that  scattering  responsible  fcr  coda  formation 
is  a  relatively  uniformly  random  process  ^^ith  only  minor 
subtle  differences  between  the  different  locations 
evaluated  in  this  study. 

Dainty  et  al.  (1987)  report  turbidity  values  for 
local  events  in  South  Carolina  that  range  in  magnitude  from 
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Figure  3.  (A/t^)  model,  turbidity  versus  frequency  by 
common  component  for  all  events  and  stations.  The 
least-square  best  fit  line  is  shown.  The  symbols  represent 
the  three  events:  Event  1  (+) ,  Event  2  (x)  and  Event  3  (o) 
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Figure  4.  (A/t^)  model,  turbidity  versus  frequency  by 
common  event  tor  all  stations  and  components.  The  least- 
sq»iare  best  fix,  line  is  shown.  Symbols  represent  the  three 
components:  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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RSON 


RSCP 


Figure  5.  (A/t^)  model,  turbidity  versus  frequency  by 
common  receiver  for  all  events  and  components.  The  least- 
square  best  fit  line  is  shown.  Symbols  represent  the  three 
components;  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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. 007  to  .1  at  3  Hz .  These  values  determined  from  local 
events  are  in  agreement  with  this  study,  where  the  values 
were  determined  at  teleseismic  distances. 

The  source  term  is  an  estimate  of  the  source  wavelet 
spectrum.  Prior  to  analysis  the  spectrum  of  the  pre-event 
ambient  noise  was  subtracted  from  the  coda  spectrum.  The 
shape  of  the  source  term  results  will  therefore  be  a 
signal-to-noise  ratio  distribution.  Shown  in  Figure  6  for 
the  (A/t^)  model  are  the  source  term  results,  spectral 
power  (db  down  from  maximum)  versus  frequency,  plotted  by 
common  component.  The  source  term  results  appear  to  be 
similar  for  each  component. 

Source  spectrum  results  plotted  by  common  event  are 
given  in  Figure  7.  The  displays  are  similar  to  those 
plotted  by  component  but  show  less  scatter  when  plotted  by 
common  event.  If  the  source  spectral  estimates  represent  a 
true  measure  of  source  properties  these  displays  should 
have  better  continuity  than  those  sorted  by  component  or 
receiving  station.  The  two  explosions,  events  1  (m]3=6.1) 
and  2  (mj3=6.1)  appear  to  have  a  faster  high  frequency  decay 
than  event  3,  the  smaller  magnitude  (m)3=5.8)  earthquake. 
These  spectral  differences  between  the  earthquake  and 
explosion  sources  are  consistent  with  other  published 
studies  (Taylor  et  al.,  1988;  Murphy  and  Bennett,  1982). 

The  faster  decay  of  high  frequency  energy  for  a 
shallow  explosion  is  attributed  to  passage  through  low  Q 
material  both  at  the  source  and  receive  locations.  Energy 
from  the  559  km  deep  earthquake  used  in  this  study  would 
only  pass  through  low  Q  crust  near  the  receiver.  This 
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Figure  6.  (A/t^)  model,  source  term  solution  by  common 
component  for  all  events  and  stations.  Values  plotted  are 
db  down  from  maximum.  Symbols  represent  the  three  events: 
Event  1  (+) ,  Event  2  (x)  and  Event  3  (o) . 
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Figure  7.  (A/t^)  model,  source  term  solution  by  common 
event  for  all  stations  and  components.  Values  plotted  are 
db  down  from  maximum.  Symbols  represent  the  three 
components:  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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spectral  characteristic  is  less  apparent  for  events  with 
magnitudes  above  4.5  to  5.0  (Taylor  et  al.,  1988).  The  fact 
that  the  effect  is  clearly  observed  on  magnitude  5.8  and 
6.1  events  in  this  study  is  quite  interesting  and  merits 
extended  investigation  with  a  larger  data  set.  One  possible 
explanation  is  that  the  source  spectrum  computed  for  500  to 
700  seconds  of  coda  is  more  stable  than  one  computed  from 
the  direct  arrival  alone.  Alternatively  the  source-to- 
receiver  distance  may  be  an  important  factor,  teleseismic 
in  this  study  versus  regional  for  Taylor  et  al.,  (1988). 

Source  depth  can  also  affect  the  decay  of  energy  at 
high  frequency  (Taylor  et  al.,  1988).  Deep  explosion  models 
involving  a  sudden  pressurization  of  a  spherical  cavity 
result  in  an  (i)”2  spectral  decay  at  high  frequencies  (cf. 
Sharpe,  1942  and  Aki  ana  Richards,  1980) .  Shallov/ 
explosions  involving  nonlinear  effects  such  as  crushing  and 
plastic  flow  can  result  in  an  «•*“ "  decay  at  high 
frequencies.  Superimposed  on  figure  7  for  the  shallow 
events  1  and  2  is  the  uT^  model  and  for  the  deep  event  3 
the  model.  These  models  appear  to  fit  i.hc  high 
frequency  decay  of  energy  observed  in  the  solutions.  Due  to 
the  high  amount  of  scatter  and  the  limited  data  set 
analyzed  it  is  difficult  to  draw  any  definite  conclusions 
from  this  comparison.  This  comparison  would  be  more 
meaningful  in  a  study  with  a  larger  data  set. 

For  the  sake  of  completeness  Figure  8  shows  the 
source  spectrum  results  from  the  three  events  plotted  by 
common  receiving  station.  The  motivation  for  this  sorting 
of  the  source  term  results  was  to  determine  if  near- 
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Figure  8.  (A/t^)  model,  source  term  solution  by  common 
receiver  for  all  events  and  components.  Values  plotted  are 
db  down  from  maximum.  Symbols  represent  the  three 
components:  vertical  (+) ,  radial  (x)  and  transverse  (o) . 
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receiver  effects  were  affecting  the  source  spectral 
estimates  and  therefore  corrupting  the  source  information. 
No  strong  receiver  dependence  is  observed  for  the  source 
term. 

(A/t+B)exp(-Ct)  MODEL  RESULTS 

It  is  not  possible  to  separate  the  source  terra  from 
turbidity  with  the  (A/t)  model  results.  In  this  case  the 
product  of  turbidity  and  source  spectra  (TS)  is  shown  in 
the  displays  that  follow.  The  same  average  crustal  velocity 
used  in  the  analysis  of  the  (A/t^)  results  is  used  here 
(3.5  km/s) .  The  value  of  crustal  thickness  used  in  the 
calculations  is  35.0  km.  The  TS  results  shown  below  were 
determined  from  the  sum  of  the  A  and  B  model  parameters. 
When  TS  was  computed  from  A  and  B  independently,  the 
results  had  similar  shapes  but  magnitude  differed  by  as 
much  as  1000.  These  magnitude  differences  suggest  that 
turbidity  may  be  different  for  the  two  types  for  scattering 
represented  by  the  plane  wave  model. 

Shown  in  Figure  9  is  log(TS)  versus  frequency  by 
common  component.  The  displays  are  similar  in  each  case 
suggesting  that  each  component  is  measuring  similar  earth 
and  source  information.  The  least-square  best  fit  line  is 
also  shown  in  TS  result  displays.  Although  the  shape 
suggests  that  the  source  term  is  influencing  the  results 
and  thus  a  best  fit  line  may  not  be  meaningful,  the  lines 
were  left  because  they  act  as  reference  points  for 
comparisons . 

Figure  10  is  log(TS)  versus  frequency  plotted  by 
common  event.  Large  variances  are  still  observed  in  the 
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Figure  9.  (A/t)  model,  product  of  turbidity  and  source  term 
versus  frequency  by  common  component.  Symbols  represent  the 
three  events,  event  1  (+) ,  event  2  (X)  and  event  3  (0) . 
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Figure  10.  (A/t)  model,  turbidity  and  source  term  product 
versus  frequency  by  common  event.  Symbols  represent  the 
3-components,  vertical  (+) ,  radial  (X)  and  transverse  (0) . 
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solutions  with  this  sorting.  It's  Interesting  that  the 
results  for  the  two  explosions  (event  1  and  2)  have  strong 
similarities.  The  decay  of  energy  at  high  frequency  appears 
to  be  slower  for  the  earthquake  than  the  explosions.  If 
turbidity  is  relatively  constant  with  frequency  these 
shapes  could  represent  estimates  of  the  source  spectra. 

Figure  ll  is  log(TS)  versus  frequency  plotted  by 
common  receiving  station.  There  appears  to  be  less 
variability  in  the  solutions  when  sorted  by  common 
receiving  station.  This  again  suggests  that  processes  near 
the  receiving  station  are  dominant  contributors  to  these 
solutions.  While  this  is  reasonable  for  the  turbidity  term 
it  is  not  expected  for  the  source  term. 

An  additional  comparison  between  the  two  models  was 
made  by  computing  the  product  of  turbidity  and  source  term 
(TS)  for  the  (A/t^)  model.  In  Figure  12  the  (A/t^)  and 
(A/t)  model  TS  solutions  are  plotted.  In  all  cases  the 
results  are  displayed  as  db  down  from  maximum  versus 
frequency.  The  solutions  are  for  each  of  the  three  events 
at  station  RSON.  The  (A/t)  model  solutions  appear  to  have 
slightly  lower  variability  than  the  (A/t^)  model  results. 
There  is  however  good  agreement  between  the  two  sets  of 
solutions.  For  example,  the  corner  frequency  for  each  event 
occurs  at  approximately  the  same  frequency  (1.5-2. 0  Hz)  for 
both  models.  The  fall  off  of  energy  at  high  frequency  is 
also  similar.  For  each  event  a  reference  line,  which 
follows  the  high  frequency  decay  ot  energy,  is  drawn  for 
easier  comparison.  These  similarities  suggest  that  the 
two-term  coda  analysis  approach  is  not  as  model  sensitive 
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Figure  ll.  (A/t)  model,  product  of  turbidity  and  source 
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Figure  12.  Product  of  turbidity  and  source  term  versus 
frequency  for  both  the  (A/t^)  and  (A/t)  models.  Comparison 
of  all  three  events  at  station  RSON,  Symbols  represent  the 
3-ccmponents,  vertical  (+) ,  radial  (x)  and  transverse  (o) . 


199 


as  was  originally  thought 


CONCLUSIONS 

The  results  of  this  feasibility  study  are 
encouraging.  The  use  of  two-term  models  in  coda  analysis 
was  shown  to  be  a  valid  approach  to  obtaining  source  and 
earth  property  estimates.  Source  spectrum  and  earth 
turbidity  information  was  obtained  from  the  A  and  B  model 
parameters . 

Turbidity  estimates  were  found  to  be  in  agreement 
with  other  estimates  in  the  literature  (Dainty  et  al., 

1987)  and  ranged  in  value  from  .007  to  0.2.  Turbidity  in 
this  study  appears  to  be  dominated  by  near-receiver 
properties.  Much  lower  variability  was  observed  in  the 
results  when  they  were  sorted  by  common  receiver  than  by 
common  component  or  event.  Turbidity  estimates  were  only 
possible  with  the  (A/t^)  model.  They  were  found  to  be 
nearly  constant  with  frequency. 

Although  turbidity  could  not  be  uniquely  separated 
with  the  (A/t)  model,  the  source  spectrum  shape  was 
estimated  by  assuming  that  turbidity  was  constant  with 
frequency.  Both  the  (A/t^)  and  (A/t)  model  results  show 
that  all  three  components  (vertical,  radial  and  transverse) 
contain  similar  information.  The  explosion  events  had 
faster  decay  of  energy  at  high  frequency  than  the 
earthquake.  Source  spectrum  estimation  from  coda  analysis 
was  shown  to  be  a  potentially  new  way  of  obtaining  spectral 
information  about  the  source  that  could  be  useful  in 
discriminating  between  earthquakes  and  explosions. 

The  (A/t)  model  solutions  were  found  to  have  a  strong 
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near-receiver  component.  Better  clustering  and  lower 
variability  was  observed  when  the  (A/t)  model  results  were 
plotted  by  common  receiver.  When  the  solutions  were  sorted 
by  common  event  the  explosion  results  were  observed  to 
decay  faster  at  high  frequency  than  the  earthquake.  This 
type  of  difference  has  been  reported  in  the  literature  as  a 
means  of  discriminating  between  explosions  and  earthquakes 
at  local  and  regional  distances. 

In  paper  1  lower  variability  was  observed  in  the 
coda-Q  solutions  with  the  (A/t^)  model  results  suggesting 
that  this  model  fit  the  data  better.  In  this  paper  the 
(A/t)  model  solutions  were  observed  to  have  the  same  or 
slightly  lower  variability  than  the  (A/t^)  model.  Based  on 
the  results  of  this  study,  including  papers  1  and  2,  it  is 
not  possible  to  conclude  that  the  (A/t^)  or  (A/t)  model  was 
best.  On  the  contrary  the  results  suggest  that  the  two-term 
coda  analysis  approach  is  not  as  model  sensitive  as  was 
originally  suspected  and  either  model  will  generate  similar 
solutions. 

The  results  of  this  study  suggest  that  multiple 
scattering  effects  can  be  modeled  in  coda  and  additional 
earth  and  source  information  can  be  extracted  from  a  coda 
analysis.  While  these  preliminary  results  using  two-term 
models  demonstrate  the  potential  benefits  of  such  an 
approach,  they  also  suggest  that  additional  studies  with 
larger  data  sets  are  required  to  better  understand  the 
significance  of  the  solutions. 
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APPENDIX 


CODA  DUE  TO  SCATTERING  OF  A  PLANE  WAVE 
WITHIN  A  LAYER 

(By:  A.  M.  Dainty,  Earth  Resources  Laboratory,  M.I.T.) 

MODEL,  ASSUMPTIONS  AND  SIMPLIFICATIONS 
The  physical  model  is  shown  in  Figure  A-la  and 
consists  of  a  layer  of  scatterers  of  thickness  h  underneath 
the  receiver.  Where  a  more  specific  model  is  required,  the 
layer  will  be  assumed  to  be  the  crust.  Two  types  of 
propagation  will  be  considered,  propagation  in  a  half  space 
("body  wave  scattering")  and  in  the  layer  after  the  manner 
of  Lg  ("surface  wave  scattering") .  The  source  is 
considered  to  be  a  plane  wavelet,  short  in  time  but 
nonetheless  with  a  narrow  frequency  content  around  a  center 
frequency  f.  By  this  means  we  introduce  the  frequency- 
dependent  amplitude  of  the  primary  wave  So(f).  The  source 
spectrum  S  =  Sq^. 

The  following  assumptions  and  simplifications  are 

made: 

1.  The  scatterers  are  parametrized  by  a  turbidity  T(f) 
where  T(f)/(47r)  is  the  scattering  cross-section  per  unit 
volume  per  unit  solid  angle,  and  the  scattering  is 
considered  to  be  isotropic.  Then  the  scattered  amplitude 
sj^  at  distance  r^  from  a  scattering  volume  dV  is 


Figure  A-la.  Scattering  geometry,  single  scattering. 


Figure  A-lb.  Integration  geometry,  single  scattering. 
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for  body  wave  scattering  and 


So2  dV  T 

— - -  (2) 

4  jr  h  r^ 

for  surface  wave  scattering  (i.e.,  dispersion  is  ignored, 
as  appears  reasonable  for  Lg;  Herrmann  and  Kijko,  1983) . 

The  factor  of  1/h  corrects  for  the  spreading  of  the  energy 
over  the  thickness  of  the  layer. 

2.  The  total  scattered  amplitude  between  lag  times  t  and 
t+dt  after  first  arrival  time  is  given  by 


S 


(3) 


=  P(f,t)  dt 


(4) 


where  the  volume  of  integration  is  taken  between  the 
surfaces  of  constant  lag  time  t  and  t+dt  and  it  is  assumed 
that  the  relative  phase  between  scattered  arrivals  is 
random,  thus  the  squared  amplitudes  are  summed 
(integrated) .  The  relation  (4)  between  amplitude  and  coda 
power  P(f,t)  is  taken  from  Aki  and  Chouet  (1985) . 

3.  The  effect  of  attenuation  after  the  first  arrival  is 
given  by  exp(-27rft/Q)  . 

4.  The  acoustic  problem  with  a  mean  velocity  v  is  solved. 

5.  The  incident  wave  is  assumed  to  be  vertically  incident. 


SINGLE  SCATTERING 

BODY  WAVE  SCATTERING 

Refer  to  Figure  A-la  for  the  geometry  and  the 
definition  of  cylindrical  polar  coordinated  r,z,0,  where  0 
is  the  azimuth  angle.  The  lag  time  t  for  a  wave  scattered 
at  r,z  is 
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Using 


(h  -  z)  ri  h 

t  - -  +  _±  _ 

V  V  V 


rj^2  =  j-2  +  z2, 


=  (ri  -  z)  /  V 


(5) 


we  obtain  for  the  surface  of  constant  scattering  lag  time 

r2  2  t  z 

“2  “  “  ' 

V 

which  is  a  paraboloid  of  revolution  truncated  at  z  =  0  and 
z  =  h.  The  volume  of  integration  is  sketched  in  Figure 
A- lb.  The  volume  element:  dV  =  r  d0  dr  dz, 

but  we  may  immediately  integrate  over  <p  to  give 

dV  =  2  TT  r  dr  dz. 

The  differential  dr  may  be  expressed  in  terms  of  dt  by 
noting  that  at  constant  z,  from  (5)  and  Figure  A-la, 


dri  ri  ri 

dr  =  -  =  —  dr 3^  =  —  v  dt 

sin  0  r  r 


(6) 


Then  combining  (3) ,  (l)  and  (6), 


=  V  dt 


h  s^2  T 


dz  =  V  dt 


2  ri 


'h  Sq2  T 
0  2(vt+z) 


dz, 


using  (5) .  Integrating,  using  (4)  and  adding  attenuation, 


1^  h  -  2  n  f  t 

Pl(f,t)  =  -  -Sq^  T(f)  V  ln(  1  +  —  )  exp( - ) 

2  vt  Q 


(V) 


There  are  two  extreme  cases.  For  h  >  vt,  the 
constant  1  in  the  logarithm  may  be  neglected  ("thick 
layer") .  More  interesting  is  the  case  h  <  vt  ("thin 
layer") .  Then  the  logarithm  in  (7)  may  be  expanded  to  give 


_  1  -  2  TT  f  t 

Pl(f,t)  =  T  h - exp( - ) 

2  t  Q 


(8) 


which  is  identical  with  Equation  (10)  of  Dainty  (1985) 
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allowing  for  a  difference  in  the  definition  of  turbidity. 

In  fact,  this  formula  may  be  derived  directly  by  means  of 
the  approximations  r^^r,  vt=r  and  dV=hdA=2jrhrdr=27rhrvdt 
used  in  (1)  and  (3).  Physically,  this  is  equivalent  to 
saying  that  the  depth  of  the  scatterer  may  be  neglected  and 
the  scattering  volume  approximated  by  the  volume  between 
two  cylinders  of  depth  h;  dA  is  the  area  on  the  surface 
between  the  two  cylinders.  This  will  be  referred  to  as  the 
"thin  layer  approximation”.  Figure  A-lc  compares  h/vt  with 
ln(l+h/vt);  we  see  that  for  vt/h  >  1,  the  thin  layer 
approximation  is  quite  accurate. 

SURFACE  WAVE  SCATTERING 

First,  we  note  for  Lg  scattering  the  thin  layer 
approximation  must  hold  because  the  scatterer  must  be 
sufficiently  far  away  to  allow  the  phase  to  develop.  Since 
Ig  is  a  collection  of  higher  modes  that  may  be  viewed  as 
po.'st-critical  crustal  S  bounces,  this  implies  distances  of 
at  least  two  crustal  thicknesses,  i.e.,  2h.  For  the 
fundamental  mode  (Rg)  there  is  no  such  limitation  on 
distance,  but  this  mode  is  confined  to  the  upper  5  km  for 
frequencies  greater  than  0.5  Hz  and  only  scatterers  in  this 
depth  range  are  important.  Subsequent  discussion  shall 
refer  to  Lg. 

Under  the  thin  layer  approximation,  we  substitute 
ri  =  r  and  dV  =  27rhrdr  =  25rhrvdt  in  (2)  and  (3)  and  using 
(4)  obtain: 

1  -  -  2  ^  f  t 

=  -  Sq^  t  V  exp( - )  (9) 

2  Q 

This  is  identical  to  Equation  (8)  of  Dainty  (1985)  after 
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correcting  for  a  factor  of  1/h  as  discussed  in  (2) . 


DOUBLE  SCATTERING 


BODY  WAVE  SCATTERING,  THIN  LAYER  APPROXIMATION 

The  geometry  is  shown  in  Figure  A“2a  in  map  view; 
r2  +  r3  =  vt  =  2a 

whei-e  a  is  the  semimajor  axis  of  the  ellipse,  and  r^  =  2c, 
where  c  is  half  the  distance  between  foci  of  the  ellipse; 
the  receiver  is  one  focus  and  the  first  scattering  volume 
at  dA^  is  the  other.  Then  the  scattered  amplitude  at  dA2 
from  dA^  is 


2  = 


Sq^  T  h  d  hi 


4  n 

and  the  scattered  amplitude  at  the  receiver  from  double 


Sq2  t2  h2  (5^2 


scattering  is 

2  Si^  T  h  d  A2  So2  t2  h2  dAi  dA2 

^  An  16  7r2  r2^  r32 

from  (1)  and  the  thin  layer  approximation.  Thus 


r  dAi  dA2 

Sn  =  So2  t2  h2  - - - 

J  16  w2  r2^  r2^ 

The  differential  dA^  =  27rrj^dr]_.  Following  Gao  et  al. 
(1983),  the  following  relations  may  be  use: 


dA2  =  r3  dr3  de 


a2  +  c2  -  2  a  c  COS0 


a2  -  c2 


r2  =  vt-r3  - - »  ^3 

a  -  c  COS0 


a  -  c  COS0 


a2+c2-2ac  COS0  a2+c2-2ac  cos0  v  dt 

dr3  - - -  da  =  - - - 

(a-c  cos0)^  (a-c  cos0)^  2 

Using  these  relations  and  (4) , 

i  -2  n  f  t 

P2{f,t)  - - So^  t2  h2  K2  exp{ - ), 

8  n  Q 
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where 


8  JT  u  du  4  jr 

=  -  - -  =  - t  -ln(  l“Uro2  )  ] 

t  j  0  1  -  u2  t 

We  must  now  consider  the  question  of  the  proper  value 
of  Ujj.  At  first  sight  it  is  Uj^  =  1,  corresponding  to  ri  = 
vt.  This  however  leads  to  a  logarithmic  singularity  in  K2. 
This  problem  occurs  because  the  limit  r^  =  vt  is  the  case 
where  the  scattering  angle  a  (Figure  A-2a)  is  0,  tbs 
ellipse  degenerates  to  a  line,  and  there  are  second 
scatterers  infinitely  close  to  the  receiver,  leading  to  a 
singularity  of  order  in  the  integrand  of  (10) . 

However,  Sato  (1982)  points  out  that  waves  scattered  at 
angles  a  <  30°  do  not  form  separate  scattered  waves  in  the 
coda  but  interfere  with  the  incident  wave  to  cause 
fluctuations.  To  exclude  such  waves,  take  a  =  30°,  and  r2  = 
r3  to  give  rjj^  =  vt  cosl5°,Uj|j=cosl5°  =  0.965,  and  K2  = 
10.87r/t.  Then 

non!  -2;rft 

P2(f,t)  =  1.35  So2  t2  h2 - exp( - )  =  1.35  T  h  Pi 

2  t  Q 

Dainty  (1985)  found  (Th)  to  be  of  the  order  of  10“^ 
(correcting  for  the  error  in  Equation  (8))  for 
Semipalatinsk  explosions,  giving  P2  ~  0.1  P^. 

The  case  of  triple  scattering  in  this  approximation 
has  also  been  solved,  giving  (refer  to  Figure  A-2b  for  the 
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geometry) 


1  ^  -2  TT  f  t 

P  - -  So2  t2  h-i  K3  exp( - ), 

32  Jr2  Q 


where 


ir'^m  r'^m  dwf’^  d/i  f’''  1-w-B  cose  de 

K3=-  udu  —  - - -  - 

tj  0  Jo  J  0  1“u2-2w-2UW  cos  M  Jo  1+U2-2w(W+1) +2C  CC 

with  b2  =  u2+w2+2uw  cos  fi,  C  =  uw+(l-w)B  and  Ujjj=Wjjj= 0.965. 


SURFACE  WAVE  SCATTERING,  THIN  LAYER  APPROXIMATION 

Refer  again  to  Figure  A“2a  for  the  geometry.  The 
scattered  amplitude  at  dA2  from  dA^^  is 

So2  T  dAi 

Sj2  =  - - 

4  w  r2 

and  the  scattered  amplitude  at  the  receiver  from  double 
scattering  is 

2  _  Si2  T  dAj  _  So2  t2  dAi  dA2 
^  4  r3  16  7r2  r2  r3 

from  (2) .  Thus 

S2  =  3^2  t2  - - -  (11) 

J  16  T2  r3 

Expressing  dA^,  dA2,  r2  and  r3  in  the  same  manner  as  before 
and  using  (4), 

1  .,0  -2  n  f  t 

P2(f,t)  =  So2  t2  V  K2  exp( - ), 

8  TT  Q 


where 


1  r^m  r2’''  de 


K2  =  -  dr]^  - . 

2  Jo  Jo  cose 

Writing  as  before  rj^  =  uvt,  a  =  vt/2,c  =  uvt/2, 

‘Um  r2ir  de 

K2  =  vt  u  du  - . 

Jo  Jo  1 
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The  definite  integral  over  0  has  the  value  (CRC  Math 
Tables,  Integral  #  645)  2jr/y(i-u2) ,  leading  to 
K2  =  27rvt(l  -  /Tl=^)  ) . 

We  see  that  in  this  case  we  may  take  rj^  =  vt,  Uj^  =  1; 
this  is  because  the  singularity  in  (11)  is  only  of  the 
order  of  l/r2.  Then  K2  =  2  tt  v  t  and 

-2  7rft  1 

P2(f»t)  =  -  7*2  v2  t  exp( - )  =  -  T  V  t  Pi 

4  Q  2 


BODY  WAVE  SCATTERING,  THICK  LAYER  APPROXIMATION 
Refer  to  Figure  A-3  for  the  geometry. 

1  00  -2  TT  f  t 

P2(f»t)  - - t2  V  K2  exp( - ) 

Bn  Q 


where 


K2  =  2jrvt 


•u 


m 

udu 


'’*^/2  sin  0  d0  (1+u  cos0)2+u2  sin0 
- ln( - ] 

0  1+u  cos  0  (1+u  cos0)2-u2  sin0 


REMARKS 

It  will  be  seen  that  in  the  thin  layer  approximation, 
for  body  wave  scattering 


A  “2  jr  f  t 

Pb(f/t)  =  -  exp( - ) 

t  Q 

for  single  scattering  and  all  higher  orders  of  scattering. 
For  surface  wave  scattering. 


-2  jr  f  t 

Ps(f»t)  =  (  B  +  C  t  +  •••  )  exp(  -  ) 

Q 

where  B  is  single  scattering,  C  is  double  scattering. 


213 


Figure  A-3.  Scattering  geometry,  double  scattering,  thick 
layer. 
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